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Abstract 

Hydrocodes  are  valuable  tools  in  the  modeling  of  shock  wave  propagation  through 
solids  due  to  high  speed  impact  phenomena.  To  model  high  speed  impacts,  hy¬ 
drocodes  use  the  conservation  equations  of  mass,  momentum  and  energy,  constitutive 
equations  to  model  the  materials  response  and  equations  of  state  (EOS).  CTH  is  a 
hydrocode  built  with  the  ability  to  use  multiple  EOSs,  including  the  semi-emperical 
Mie-Gruneisen  EOS  and  tabular  Sesame  EOS.  Modeling  high  speed  impacts  necessi¬ 
tates  modeling  the  non-equilibrium  thermodynamic  states  caused  by  these  impacts. 
A  discussion  of  the  non-equilibrium  thermodynamics  that  may  be  applied  to  the  re¬ 
gion  directly  behind  a  shock  wave  is  presented,  including  details  of  recent  attempts  to 
model  non-equilibrium  impact  phenomena  in  solids.  Also,  in  order  to  better  determine 
the  applicability  of  the  Mie-Gruneisen  EOS  and  the  Sesame  EOS  in  situations  that  in¬ 
clude  non-equilibrium  thermodynamics,  the  high  speed,  uniaxial  impacts  between  two 
iron  bars  are  modeled  in  CTH.  These  impacts  are  modeled  using  the  Mie-Gruneisen 
EOS,  Sesame  EOS  and  a  two  state  EOS  called  PTRAN.  The  results  generated  using 
these  different  EOS  are  then  compared  to  experimental  data  to  determine  how  well  the 
different  EOS  model  the  thermodynamic  non-equilibrium  physics  of  the  high  speed 
impacts.  The  differences  between  the  Mie-Gruneisen  EOS  and  the  Sesame  EOS  are 
established.  A  finite  volume  uniaxial  hydrocode  is  validated.  Finally,  CTH  is  shown 
to  be  able  to  model  some  irreversibilities  occurring  in  impact  phenomena. 
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Comparison  of  Thermodynamic 


Equilibrium  and  Non-Equilibrium 
Representations  of  Materials 

I.  Introduction 

1 . 1  Introduction 

A  hydrocode  is  a  computer  code  that  is  able  to  model  the  waves  caused  by 
high  speed  impacts  as  they  travel  through  solid  materials.  As  Zukas  [43]  explains 
in  his  hydrocode  book,  early  hydrocodes  were  given  the  name  because  they  relied 
on  the  hydrodynamic  assumption.  The  hydrodynamic  assumption  states  that  if  the 
pressures  generated  in  a  high  speed  impact  are  much  greater  than  the  strength  of  the 
materials  involved  in  the  impact,  then  the  strength  may  be  ignored  and  the  materials 
may  be  modeled  as  fluids.  Modern  hydrocodes  are  able  to  account  for  the  strength  of 
the  materials  by  using  constitutive  equations  in  these  impact  scenarios,  however  the 
name  hydrocode  has  stuck. 

Today,  hydrocodes  are  used  to  model  an  array  of  different  hypervelocity  impact 
scenarios.  These  scenarios  may  include  debris  striking  a  space  vessel,  or  flying  debris 
from  an  explosion  striking  a  hardened  target.  Cinnamon  [9],  Szmerekovsky  [35]  and 
Laird  [24]  used  the  Sandia  National  Labortories  hydrocode  CTH  to  model  the  gouging 
caused  by  impact  between  the  shoes  of  the  sleds  used  at  the  Holloman  High  Speed 
Test  Track  (HHSTT)  and  the  railing  of  the  track  itself. 

Hydrocodes  use  three  different  types  of  equations  to  model  impacts.  These 
time  dependent  equations  include  conservation  equations,  the  previously  mentioned 
constitutive  equations  and  equations  of  state  [43]. 

According  to  continuum  theory,  the  conservation  equations  include  the  con¬ 
servation  of  mass,  momentum  and  energy.  These  equations  give  us  the  variables  of 
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velocity  (v),  density  (p),  pressure  (P),  and  internal  energy  (e).  Temperature  (T),  may 
also  be  determined  using  these  equations,  but  is  not  necessary.  In  isothermal  sys¬ 
tems  the  temperature  and  internal  energy  of  the  system  remain  constant.  Then  the 
conservation  equations  and  an  isothermal  relationship  between  two  thermodynamic 
variables,  excluding  temperature  and  internal  energy  may  act  as  a  complete  system 
of  equations  and  thus  are  the  only  equations  necessary  to  solve  for  desired  material 
properties  [30]. 

The  previously  mentioned  impact  conditions  are  not  isothermal,  requiring  an 
additional  equation  to  close  the  system  of  equations,  hence  the  need  for  an  equation 
of  state  (EOS).  EOS  act  as  a  bridge  between  the  macroscopic  continuum  mechanics 
represented  by  the  conservation  equations  and  the  microscopic  thermodynamics  of  the 
system  [43].  Also,  the  extremely  high  pressures  generated  in  the  aforementioned  im¬ 
pact  scenarios  cause  the  EOS  used  to  dominate  the  answers  generated  by  hydrocodes. 
This  necessitates  that  the  EOS  accurately  model  the  thermodynamics. 

Unlike  in  the  study  of  gas  dynamics,  where  a  general  EOS  can  be  applied  to  a 
wide  range  of  gasses,  the  EOS  of  solids  tend  to  be  specific  to  a  material  or  material  type 
(metals,  ceramics,  etc...).  Also,  depending  on  the  origin  of  the  EOS  (empirical,  semi- 
empirical,  theoretical,  etc...),  said  EOS  may  have  a  limited  range  of  applicability  [43]. 
For  example,  empirical  based  equations  of  state  may  mean  nothing  when  used  in 
conditions  far  from  the  experiments  used  to  generate  them.  The  applicability  of 
a  theoretical  EOS  will  depend  on  the  thermodynamics  used  to  derive  it.  In  the 
area  immediately  behind  a  shock  wave,  for  example,  conditions  are  changing  at  an 
extremely  high  rate,  so  there  exists  a  state  of  non-equilibrium  necessitating  the  use  of 
non-equilibrium  thermodynamics  (NET).  Conversely,  in  regions  where  the  conditions 
are  not  changing  as  fast,  a  state  of  local  equilibrium  may  be  assumed  allowing  for  the 
use  of  equilibrium  thermodynamics  [26]. 
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Equilibrium  based  thermodynamics  were  used  in  the  generation  of  the  popular 
Mie-Gruneisen  EOS  and  in  the  generation  of  the  tabular  EOS  such  as  the  Sesame 
EOS  [22], 

In  the  1930’s  work  began  in  Classical  Irreversible  Thermodynamics  (CIT),  which 
used  thermodynamic  fluxes  in  addition  to  the  classic  thermodynamic  state  variables 
in  the  description  of  a  system  in  an  assumed  local  equilibrium  state.  The  assumption 
of  a  local  equilibrium  state  allows  CIT  to  extend  beyond  the  range  of  equilibrium 
thermodynamics,  but  still  requires  a  system  to  remain  ’’near”  an  equilibrium  state  [26]. 

In  the  1960’s  another  form  of  irreversible  thermodynamics  referred  to  as  ratio¬ 
nal  dynamics  (RT)  was  introduced.  RT  assumes  materials  have  memories  and  that 
complimentary  variables  (entropy,  internal  energy,  heat  flux,  etc...)  are  functions  of 
the  histories  of  the  independent  state  variables  [27] . 

Two  branches  of  modern  irreversible  thermodynamics  are  referred  to  as  ex¬ 
tended  irreversible  thermodynamics  (EIT)  and  internal  variable  theory  (IVT).  EIT 
uses  macroscopic  thermodynamic  fluxes  as  non-equilibrium  state  variables  to  com¬ 
plement  the  use  of  classical  equilibrium  state  variables  in  the  description  of  a  local 
non-equilibrium  state  [18].  IVT  identifies  additional  variables  that  describe  the  mi¬ 
croscopic  effects  of  a  system. 

Today  researchers  such  as  Lu  and  Hanagud  [26]  [27]  are  able  to  use  a  mixture  of 
irreversible  thermodynamic  theories  to  model  the  non-equilibrium  thermodynamics 
of  a  solid  material’s  response  to  impact  loading. 

1 . 2  Procedure 

There  are  several  portions  to  the  numerical  investigation  in  this  investigation. 
These  parts  include  the  comparison  of  the  results  generated  using  the  Mie-Gruneisen 
and  PTRAN  EOS  to  those  produced  using  the  Sesame  EOS,  the  validation  of  the  finite 
volume  Eulerian  hydrocode  and  a  comparison  of  the  results  generated  using  CTH  to 
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those  generated  using  the  Non-Equilibrium  Thermodynamic  hydrocode  developed  by 
Lu  and  Hanagud. 

Each  portion  of  the  numerical  investigation  composes  its  own  chapter  to  avoid 
confusion.  Each  chapter  consists  of  a  description  of  the  numerical  procedure  followed 
immediately  by  the  results  generated.  Chapter  2  discusses  the  theories  behind  the 
modeling  of  equilibrium  and  non-equilibrium  thermodynamic  impact  systems.  Chap¬ 
ter  3  provides  the  procedures  and  results  from  comparing  the  Mie-Gruneisen  and 
PTRAN  EOS  to  the  Sesame  EOS  using  uniaxial  impacts  of  aluminum  and  iron,  a 
reversibility  test  to  determine  if  the  Sesame  EOS  can  model  irreversibilites  in  the 
regions  of  a  phase  change,  and  a  recreation  of  the  flyer  plate  experiments  performed 
by  Cinnamon  [9]  using  1080  steel.  Chapter  4  explains  how  the  finite  volume  hy¬ 
drocode  is  validated  comparing  the  results  generated  by  the  finite  volume  hydrocode 
to  those  generated  by  CTH.  Chapter  5  compares  the  results  generated  by  CTH  to 
those  produced  by  Lu  and  Hanagud’s  non-equilibrium  hydrocode.  Finally,  Chapter  6 
presents  any  conclusions  made  and  provides  recommendations  for  the  further  study  of 
non-equilibrium  thermodynamics  and  the  development  of  the  finite  volume  Eulerian 
hydrocode. 
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II.  Theory 

This  chapter  discusses  the  basic  concepts  behind  the  numerical  modeling  of  dynamic 
impact  phenomena.  A  general  overview  of  the  numerics  behind  the  hydrocodes  is 
provided  along  with  the  basic  form  of  the  equations  used  in  these  hydrocodes.  A 
Thermodynamics  review  is  then  presented  followed  by  the  development  of  equations 
of  state.  This  will  be  followed  by  an  introduction  to  current  methods  used  to  model 
the  non-equilibrium  thermodynamics  in  impact  scenarios.  All  the  numeric  simulations 
in  this  investigation  are  performed  using  a  uniaxial  strain  solver  in  CTH,  therefore 
all  the  equations  in  this  theory  section  are  written  in  uniaxial  form  where  applicable. 

2. 1  Numerical  Modeling 

Two  programs,  the  hydrocode  CTH  and  a  finite  volume  code  written  in  Fortran 
90,  are  used  to  compare  the  abilities  of  different  equations  of  state  to  model  high  speed 
impact  phenomena.  A  brief  description  of  how  each  program  works  is  presented  here 
to  familiarize  the  reader  with  these  two  programs. 

2.1.1  Numeric  Terminology.  Two  terms  used  while  analyzing  a  code’s  abil¬ 
ity  to  numerically  model  discontinuities  such  as  a  shock  wave  (defined  in  Section 
2.3)  are  dissipation  and  dispersion.  Dissipation  refers  to  the  process  of  taking  a  dis¬ 
continuity  and  spreading  it  out  over  several  numerical  grid  cell  widths.  Dispersion 
refers  to  oscillations  in  the  numerical  solution  caused  by  discontinuities,  and  is  often 
located  around  the  discontinuity  itself.  Figure  2.1  illustrates  the  difference  between 
the  types  of  numerical  errors  compared  to  an  ’’ideal”  representation  of  a  numerical 
discontinuity  [36]. 


a.  Dissipation  b.  Discontinuity  c.  Dispersion 

Figure  2.1:  Illustration  of  Dissipation  and  Dispersion  Error 


5 


2.1.2  CTH.  Hertel  et  al.  [15]  provide  a  concise  explanation  of  how  CTH 
works.  CTH  is  an  Eulerian  based  hydrocode  that  uses  a  two  step  process  to  solve  the 
conservation  equations.  The  first  step  is  the  Lagrangian  step,  in  which  the  mesh  is 
deformed  and  the  governing  equations  are  integrated  through  time.  Mass,  momentum, 
energy  and  volume  are  conserved  during  the  Lagrangian  step.  Mass  is  conserved,  due 
to  the  movement  of  the  mesh,  which  stops  all  matter  from  moving  across  the  cell 
boundaries.  The  rest  of  the  conservation  equations  are  replaced  by  explicit  finite 
volume  forms  of  the  original  integral  equations. 


The  Lagrangian  step  is  followed  by  a  remap  step  in  which,  the  distorted  mesh  is 
re-mapped  onto  the  original  Eulerian  mesh.  During  the  remap  step  mass,  momentum, 
energy  and  volume  are  transported  from  the  deformed  mesh  to  the  original  mesh  using 
advection  algorithms  [6].  To  do  this,  CTH  first  calculates  the  flux  of  the  volume 
between  the  old  and  new  cells.  Then,  an  interface  tracking  algorithm  decides  which 
materials  from  the  old  cells  get  moved  with  the  volume  flux.  The  mass  and  energy 
of  materials  being  moved  are  then  transferred  from  the  old  cells  to  the  new  ones. 
Finally  information  from  the  interface  tracker  is  used  to  move  the  momentum  and 
kinetic  energy  of  the  cells  [15]. 


Benson  [6]  provides  a  good  illustration  of  what  is  happening  numerically  when 
the  operator  split  is  being  used.  When  modeling  a  continuum,  the  conservation  equa¬ 
tions  equations  take  the  form: 


d(J)  d(j) 

m+cai  =  l 


(2.1) 


This  equation  is  referred  to  as  the  linear  advection  equation,  where  4>  is  the  field 
variable,  c  is  the  constant  flow  velocity  and  /  is  a  source  term.  The  operator  split 
divides  this  equation  into  two  parts.  The  Langrangian  equation,  which  is  solved  by 
moving  the  mesh  is: 


d& 

m 


(2.2) 
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and  the  Eulerian  equation  solved  during  the  remap  step  has  the  form: 


dcf)  dcj) 
~dt+C~di 


=  0 


(2.3) 


In  order  to  control  the  discontinuity  caused  by  shock  waves  in  the  Lagrangian 
time  step,  a  three  term  ’’artificial  viscosity”  is  utilized.  Artificial  viscosities  are  numer¬ 
ical  tools,  added  to  the  pressures,  that  smear  the  discontinuities  in  continuum  based 
codes  over  several  mesh  widths,  allowing  hydrocodes  to  handle  any  discontinuities 
that  may  arise.  Though  the  artificial  viscosity  distorts  the  area  around  discontinu¬ 
ities,  the  solution  is  only  affected  in  the  area  of  the  shock  front  and  the  accuracy  of 
the  calculations  are  preserved.  [43] 

In  CTH,  a  vector  subset  of  the  full  viscosity  tensor  including  linear  and  quadratic 
terms  is  used.  In  three  dimensions,  the  vector  includes  xx,  yy  and  zz  terms.  A  third 
linear  term  controling  a  singular  point  during  the  update  of  the  stress  deviators  along 
the  axis  of  symmetry  when  the  two-dimensional  cylindrical  geometry  is  also  used  [15]. 


2.1.3  The  Vanderhyde  Finite  Volume  Code  (VFVC).  A  finite  volume,  Eu¬ 
lerian  hydrocode  that  utilizes  the  Mie-Gruneisen  EOS  is  also  used  to  model  uniaxial 
impacts.  This  code,  refered  to  as  the  Vanderhyde  Finite  Volume  Code  (VFVC),  was 
adopted  from  a  computational  fluid  dynamics  (CFD)  shock  tube  code.  It  solves  the 
partial  differential  equation,  shown  in  finite  difference  form,  of  [38]: 


dQ  vec  9Fvec  f(/^i  \ 

dt  +  dx  =  IWvec) 


(2.4) 


where  Qvec  is  the  vector  of  conserved  variables,  Fvec  is  the  flux  vector  and  f(Qvec )  is 
a  source  term  vector. 

For  clarification,  a  finite  volume  solver  applies  the  conservation  equations  to  a 
fixed  region  of  space,  called  a  control  volume  [36].  In  finite  volume  form  this  is  written 
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Figure  2.2:  One  Dimensional  Grid 
as: 

dVol+  f  Fvec  ■  dA  =  f  f(Qvec )  (2.5) 

J  A  JVol 

where  A  is  referred  to  as  an  area  vector  and  Vol  is  the  volume  of  the  cell,  which  in  one 
dimension  is  the  width  of  the  cell,  or  dx.  Upon  integrating  Equation  2.5  becomes: 

=  Voi  E  • dA  +  /(c,ec)  (2.6) 

The  flux  vector  numerically  represents  the  flux  of  the  conservation  variables  (mass, 
momentum,  energy,  etc)  at  the  boundary  of  the  control  volume,  and  is  used  in  the 
spatial  integration  of  the  conservation  equations.  Finally,  the  source  term  refers  to 
the  non-hyperbolic  portion  (the  right  side)  of  Equations  (2.4)and  (2.5).  An  example 
of  a  one  dimensional  grid  is  shown  in  Figure  2.2. 

The  VFVC  uses  a  flux  routine  to  solve  for  the  flux  variables  moving  through 
the  mesh  and  to  handle  any  discontinuities  that  may  present  themselves.  Fi  —  Fi+ 1 
represents  the  flux  through  the  faces  of  individual  cells  where  i  represents  the  node 
numbers  of  the  cell.  The  flux  routine  used  in  the  VFVC  is  referred  to  as  the  Lax- 
Friedrichs  flux  routine  [38].  Equation  (2.7)  gives  the  form  of  the  Lax-Friecrichs  flux. 
Lax-Friedrichs  is  a  first  order  approximation  of  the  flux  vector  values  at  the  cell  face 
and  works  by  taking  the  average  values  of  the  flux  vector  at  the  cell  boundaries  and 
adds  a  dissipation  term,  ai+i(q(xi )  —  q(xi+ 1)),  to  the  average,  allowing  the  code  to 
handle  discontinuities. 
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(2.7) 


fy(x,))  +  fy(xi+1))  +  al+i(q(Xi)  -  9(14+1)) 


A  second  order  MUSCL  Scheme  [36]  is  used  to  remove  excess  dissipation  intro¬ 
duced  by  the  Lax-Friedrich  flux  routine  and  provides  second  order  spatial  accuracy 
with  the  use  of  numerical  limiters  that  remove  non-physical  dispersion  error. 

After  the  flux  routines  solve  for  the  flux  vector  values,  the  flux  vector  values 
are  inserted  into  the  right  side  of  Equation  (2.6),  also  called  the  right  hand  vector, 
R(Qvec) ;  which  becomes: 


R(Qvec ) 


F,  -  Fi 


i+ 1 


dx 


+  D(Qvec ) 


(2.8) 


where  dx  is  the  cell  width  and  D  is  the  source  term  operator  [38] . 


Explicit  time  integration  schemes  are  used  to  integrate  the  following  PDE  with 
respect  to  time.  [38]: 


dQ 


vec 


dt 


R(Q  vec ) 


(2.9) 


Three  different  time  integration  schemes  are  used  to  integrate  through  time:  a  first 
order  Euler  Explicit  Scheme  [36] ,  a  third  order  Runge-Kutta  scheme  [38]  and  a  fourth 
order  Runge-Kutta  scheme  [36] .  The  Euler  Explicit  Scheme  used  is  given  in  Equation 
(2.10),  while  the  other  integration  schemes  are  contained  in  Appendix  1.  The  Euler 
Explicit  method  uses  first  order  forward  difference  schemes  in  space  (i)  and  time 
(n),  to  integrate  through  time.  A  more  detailed  development  of  the  uniaxial  code  is 
available  in  Appendix  1. 


Q-+l  -  Qi  ,  Q?+ 1  -  Qri 


At 


+  c- 


Ax 


=  0 


(2.10) 


2. 2  Equations  of  Motion 

In  the  numerical  modeling  of  a  dynamic  system  several  classes  of  the  equations 
of  motion  equations  are  needed.  These  equations  of  motion  include  the  conservation 


9 


equations,  an  equation  of  state  (EOS)  and  when  modeling  systems  where  the  strength 
of  materials  needs  to  be  taken  into  account,  a  set  of  constitutive  equations. 

2.2.1  Conservation  Equations.  The  conservation  equations  typically  are 
used  to  solve  for  the  values  of  the  mass,  momentum  and  energy.  Since  all  the  numer¬ 
ical  simulations  being  considered  are  uniaxial,  the  conservation  equations  below  are 
written  one  dimensionally  in  the  x  direction.  CTH  and  the  FVC  are  both  Eulerian 
based  hydrocodes,  meaning  they  track  matter  as  it  proceeds  through  a  fixed  mesh.  Eu¬ 
lerian  hydrocodes  are  typically  used  when  the  distortions  caused  by  an  event  are  very 
large,  something  pure  Langrangian  codes  have  trouble  modeling  because  cell  volumes 
may  approach  zero,  or  the  cells  become  inverted  producing  negative  volumes  [43]. 
Since  both  hydrocodes  are  Eulerian,  the  conservation  equations  are  presented  here  in 
Eulerian  form. 


2. 2. 1.1  Conservation  of  Mass.  Simply  stated  the  conservation  of  mass 
dictates  that  mass  can  neither  be  created  or  destroyed  as  a  system  undergoes  changes. 
[1] 

The  mathematical  form  of  the  conservation  of  mass  is  written  as  [37]: 


dp  d(pv) 
dt  dx 


(2.11) 


The  first  term  of  the  conservation  of  mass  accounts  for  compressibility.  Compressibil¬ 
ity  effects  refer  to  a  changes,  over  time,  in  the  density  of  a  material  at  a  given  point 
in  space,  or  a  change  of  the  entire  volume  of  a  system.  The  second  term  accounts  for 
the  change  in  the  amount  of  material  flowing  through  a  control  volume  with  respect 
to  spatial  position. 


2.2. 1.2  Conservation  of  Momentum.  Conservation  of  momentum  can 
be  stated  as:  the  time  rate  of  change  of  momentum  of  a  body  equals  the  net  force 
exerted  on  said  body.  This  is  also  Newton’s  second  law  stating  that  force  exerted  on  a 
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body  is  equal  to  the  product  of  the  mass  of  the  body  and  the  amount  of  acceleration 
the  body  is  undergoing,  or  F  =  ma  [1], 

The  Eulerian  form  of  the  conservation  of  momentum  is  written  as  [37]: 


pF  = 


d(pv)  d(pv2  —  a) 
dt  dx 


(2.12) 


where  F  is  the  force  per  unit  mass  acting  on  the  system  and  a  represents  the  stress 
being  applied  to  the  system.  In  the  conservation  of  momentum,  the  pF  term  repre¬ 
sents  the  entire  force  acting  on  the  system.  The  first  component  of  the  right  hand 
side  of  Equation  (2.12),  accounts  for  changes  in  the  amount  of  mass  flowing  through  a 
control  volume  over  time.  The  second  term  accounts  for  changes  of  the  total  momen¬ 
tum  inside  a  control  volume,  with  respect  to  spatial  position,  by  taking  the  difference 
in  net  flow  of  momentum,  pv2,  and  the  total  stress,  a,  applied  to  the  system. 


2.2. 1.3  Conservation  of  Energy.  The  conservation  of  energy  dictates 
that  energy  can  neither  be  created  or  destroyed;  it  can  only  change  form.  The  con¬ 
servation  of  energy  incorporates  the  thermodynamics  of  a  system  with  the  systems 
motion,  by  relating  a  system’s  kinetic  energy  with  other  forms  of  energy. 

The  Eulerian  form  of  the  conservation  of  energy  is  [37]: 

dp(e+\v2)  d(vp  \e  +  lv2]  —  av) 

+ =  V  +  <4 — 1  <213) 

where  q  represents  heat  diffusion.  In  English,  the  conservation  of  energy  states  the 
total  work  done,  p(q  +  Fv),  equals  the  time  rate  of  change  of  energy,  p  (e  +  |u2),  plus 
the  flow  of  energy  across  a  specified  area,  vp  [e  +  |n2] ,  minus  the  materials  resistance 
to  the  energy  flow,  av,  in  that  area  [1], 

2.2.2  Constitutive  Equations.  When  modeling  impact  dynamics,  the  re¬ 
sponses  of  the  materials  need  to  be  taken  into  consideration.  The  material  response  is 
often  represented  in  two  components,  the  hydrodynamic  response  and  the  deviatoric 
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response.  The  hydrodynamic  stress  makes  up  the  diagonal  of  the  stress  tensor,  is 
related  to  pressure  and  is  calculated  from  the  EOS  [43]. 

d hydro  ij  (2. 14) 

The  deviatoric  or  shear  stress  of  the  material  is  designated  as,  cr®.  The  total  material 
response  is  the  sum  of  the  hydrodynamic  and  deviatoric  stresses.  [38] 

aij  =  aij  ~  P5J  (2-15) 

In  high  speed  phenomena,  hydrostatic  stresses  are  typically  orders  of  magnitude 
greater  than  the  deviatoric  stress.  This  justifies  the  hydrodynamic  assumption  dis¬ 
cussed  in  the  introduction.  It  also  explains  why  the  magnitude  of  stress  produced  by 
a  very  high  speed  impact  is  approximately  equal  to  the  pressures  produced.  Through¬ 
out  this  investigation,  constitutive  relationships  are  used  to  describe  the  material  re¬ 
sponse  in  impacts.  However,  due  to  high  magnitudes  of  pressure  typically  generated, 
this  response  is  assumed  negligible. 

Two  common  constitutive  equations  used  by  CTH  are  the  Johnson-Cook  (JC) 
and  Zerilli- Armstrong  (ZA)  models.  In  this  investigation  the  availability  of  predefined 
material  coefficients  in  CTH  dictates  which  constitutive  relationship  is  used  to  model 
different  materials. 

2.2.2. 1  Johnson-Cook  Constitutive  Equations.  The  Johnson-Cook 
constitutive  model  is  an  empirical  based  relationship  that  takes  the  form: 

a  =  (A  + Ben)(l  +  Clne)(l-T*m)  (2.16) 

where  e  represents  strain,  e  represents  strain  rate,  and  A,B,n,C  and  m  are  experi¬ 
mentally  determined  material  specific  constants.  The  value  of  T*  is  found  by  the 
ratio  (T  —  Troom)  /  (Tmejt  —  Troom )  where  T  is  the  absolute  temperature  of  the  ma- 
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terial.  The  material  constants  are  determined  experimentally.  The  first  expression 
in  this  equation  represents  the  stress  as  a  function  of  strain.  The  second  and  third 
expressions  represent  strain  rate  and  temperature  effects,  respectively  [17]. 

2. 2. 2. 2  Zerilli-  Armstrong  Constitutive  Equations.  The  Zerilli- Armstrong 

constitutive  relationship  was  developed  for  the  purpose  of  improving  the  numerical 
results  generated  using  the  Johnson-Cook  relationship.  Zerilli- Armstrong  model  has 
a  strong  empirical  basis,  but  also  takes  into  account  effects  caused  by  the  microscopic 
properties  of  the  materials  such  as  grain  size  and  lattice  structure.  The  dependence 
on  the  type  of  lattice  structure  a  material  is  made  of  requires  different  forms  of  the 
Zerilli- Armstrong  relationships  [41] .  The  face  centered  cubic  (fee)  form  of  the  Zerilli- 
Armstong  relationship  is: 

<x  =  A a'G  +  c2e1/2  *  eV^T+cATine)  +  u-\  (217) 

The  body  centered  cubic  (bcc)  form  of  the  Zerilli- Armstrong  model  is: 

<X  =  A  o'G  +  Cie(-c3T+c,Tlne)  +  +  (2.18) 

For  both  the  bcc  and  fee  equations  A a'G  represents  an  additional  component  of  stress, 
k  is  the  microstructural  stress  intensity,  l is  the  inverse  of  the  square  root  of  the 
average  grain  diameter  and  clf  c2,  c3,  c4,  c5  and  n  are  experimentally  determined 
values  [41]. 


2. 2. 2. 3  Material  Values.  In  this  investigation,  the  material  response 
of  aluminum  will  be  described  using  the  Johnson-Cook  relationship,  while  iron’s  re¬ 
sponses  will  be  described  using  the  Zerilli- Armstrong  constitutive  relationship.  CTH’s 
internal  coefficients  are  used  to  describe  both  materials,  however  these  coefficients  are 
not  available.  Publicly  available  material  properties  are  given  here  for  use  in  pro- 
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grams  needing  these  input  parameters.  The  Johnson-Cook  coefficients  are  given  in 
Table  2.1  [40].  The  Zerilli-Armstrong  coefficients  are  given  in  Table  2.2  [41]. 


Table  2.1: 


Johnson-Cook  Coefficients  of  99.5%  Aluminum 


Coefficient 

Value  Units 

A 

265.0  MPa 

B 

426.0  MPa 

C 

0.015 

m 

1 

n 

0.34 

Table  2.2:  Zerilli-Armstrong  Coefficients  of  Iron 


Coefficient 

Value 

Units 

Aa'c 

0.0 

MPa 

Cl 

1033 

MPa 

c2 

na 

MPa 

c3 

0.00698 

k~l 

c4 

0.000415 

k-1 

C5 

266.0 

MPa 

n 

0.284 

k 

22.0 

MPa  mm2 

3.0 

-1 

mm  2 

2.3  Wave  Propagation  in  Solids 

This  work  depends  heavily  on  the  propagation  of  waves  through  solids,  thus  a 
brief  description  of  wave  propagation  through  solids  is  presented. 

When  an  object  is  struck  by  another,  a  compression  or  stress  wave  is  produced. 
The  type  of  stress  wave  depends  on  the  velocity  of  the  impactor  and  on  the  strength 
of  the  material  being  struck.  Three  different  wave  scenarios  are  possible.  If  the 
amount  of  strain  produced  due  to  the  applied  stress  is  below  a  material  value  called 
the  Hugoniot  Elastic  Limit  (HEL),  a  single  elastic  wave  will  propagate  through  the 
material.  If  the  amount  of  strain  produced  by  the  stress  is  above  the  HEL  an  elastic 
’’precursor”  will  precede  at  a  velocity  of  ce,  which  is  faster  than  the  wave  speed  of 
the  plastic  wave,  cp.  Figure  2.3  [43]  shows  this  second  scenario.  If  a  high  enough 
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Figure  2.3:  Stress  Wave  Profile 

value  of  stress  is  generated  by  the  impact,  the  plastic  wave  will  overtake  the  elastic 
wave  and  become  a  shockwave.  At  this  amount  of  stress,  the  material  exhibits  plastic 
behavior  and  displays  fluid  like  characteristics.  The  material  particles  behind  a  stress 
wave  move  with  a  velocity  particle  velocity,  vt  [43] . 

When  a  propagating  compression  wave  reaches  a  boundary,  how  it  reacts  is  a 
property  of  the  boundary  condition.  There  are  two  conditions  that  must  be  met  at 
the  boundary:  1)  The  forces  on  both  sides  of  the  boundary  must  be  equal,  and  2)  the 
particle  velocities  at  the  boundary  must  be  continuous.  Figure  2.4  illustrates  these 
two  conditions.  The  first  condition  leads  to  [43]: 

A\{<Ji  +  &r)  =  A2(Jt  (2-19) 

where  A  represents  the  cross-sectional  area  of  the  objects  on  both  sides  of  the  bound¬ 
ary,  and  the  subscripts  T  and  R,  represent  transmitted  and  reflected  waves,  respec¬ 
tively.  The  second  condition  states: 


Vi  —  vr  —  vt 


(2.20) 


Assuming  that  a  =  pcv  [43],  the  values  of  the  transmitted  and  reflected  stresses  can 
be  shown  to  be: 


(Jj1  —  (7  j 


2A1P2C2 

Aipici  +  A2P2C2 


(2.21) 
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Figure  2.4:  Boundary  Conditions  Between  Two  Materials 
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Figure  2.5:  Reflection  of  a  Stress  Wave 


&r  —  °7 


A2P2C2  ~  Aipici 

Aip\C\  +  A2P2C2 


(2.22) 


At  a  free  surface,  ^  approaches  0  and  p2c2  =  0  so  there  is  no  transmitted  wave 
and  gr  =  —gj.  This  means  that  a  compression  wave  reflects  as  a  tensile  or  rarefaction 
wave  (and  vice  versa),  and  that  the  particle  velocity  is  doubled  at  the  free  surface  [43]. 
Figure  2.5  shows  how  this  wave  reflection  affects  the  stress  due  to  the  compression 
wave,  and  is  included  because  free  boundaries  are  used  in  this  investigation. 

At  a  fixed  surface,  ^  approaches  infinity,  gr  =  aj  and  gt  =  0.  This  means 
that  a  compression  wave  reflects  as  a  compression  wave,  and  that  the  particle  velocity 
is  zero  at  the  free  surface  [43]. 

If  A2P2C2  =  AipiCi  then  gr  —  0  (there  is  no  reflected  wave)  and  the  value  of 
the  transmitted  stress  wave  is: 

_  F&2P2 

aT  ~  ai\j 

where  E  is  the  modulus  of  elasticity  [43] . 
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2.4  Equilibrium  and  Non- Equilibrium  Thermodynamics 

A  good  understanding  of  thermodynamics  in  necessary  to  follow  the  develop¬ 
ment  the  EOS  used  in  this  investigation.  In  this  section,  the  laws  of  thermodynamics 
are  presented.  Key  assumptions  made  in  the  use  of  equilibrium  and  non-equilibrium 
thermodynamics  are  also  presented  with  an  explanation  of  what  each  type  of  thermo¬ 
dynamics  each.  An  introduction  to  non-equilibrium  thermodynamics  is  also  presented 
here. 

2-4-1  The  Laws  of  Thermodynamics.  Anderson  [1]  simply  defines  thermody¬ 
namics  as  ’’the  science  of  energy  (and  entropy).”  The  users  of  thermodynamics  use  it 
to  approximate  the  thermodynamic  properties  of  systems,  in  other  words,  the  known 
thermodynamic  properties  of  a  system  are  used  to  calculate  its  unknown  properties 
at  equilibrium  states.  [13] 

The  thermodynamic  state  of  a  system  may  be  described  by  a  set  of  physical 
values  called  state  variables.  It  is  important  to  note  that  the  values  of  these  state 
variables  of  a  system  are  independent  of  how  the  system  arrived  at  a  particular  state. 
When  one  of  the  state  variables  changes,  the  thermodynamic  state  also  changes  [34] . 

The  first  law  of  thermodynamics  states:  ’’Energy  cannot  be  created  or  destroyed, 
only  converted  from  one  form  to  another.”  [13]  Mathematically  this  law  becomes  [39]: 


(2.24) 


de  =  dQ  —  P  *  dV 


Where  de  is  the  total  change  in  internal  energy,  dQ  is  the  amount  of  heat  the  system 
receives  from  its  surroundings  and  P  and  V  are  the  pressure  and  specific  volume  of 
the  system.  If  a  system  is  in  a  state  of  equilibrium  thermodynamics  P  and  e  are  state 
functions  that  are  related  to  the  state  variables,  specific  volume,  V,  and  absolute 
temperature  T  using  equations  of  the  form  [39]: 


P  =  P(V,  T) 


(2.25) 
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and 


e 


e(V,  T) 


(2.26) 


If  the  system  is  in  a  state  of  non-equilibrium,  additional  state  variables  that  do 
not  appear  in  the  equations  of  conservation  are  needed  to  model  the  system.  In  a 
state  of  non-equilibrium  thermodynamics,  Equations  (2.4.1)  and  (2.4.1)  become  [39]: 


P  =  P{V,  T,  N-l . Nn) 


(2.27) 


and 


e  =  e(V,  T,  Ni . Nn) 


(2.28) 


where  Nt  represent  the  additional  state  variables  needed  to  describe  the  system.  Ad¬ 
ditional  equations  are  also  needed  to  solve  for  these  variables  [39].  There  are  several 
irreversible  thermodynamic  theories  used  to  identify  and  define  these  additional  vari¬ 
ables.  Two  of  these  theories,  Internal  Variable  Theory  and  Extended  Irreversible 
Thermodynamics  are  described  in  Section  2.4.3. 

The  second  law  of  thermodynamics  deals  with  the  state  function  of  entropy. 
Gasser  and  Richards,  among  others,  define  entropy  as  ”a  measure  of  the  chaos  or 
disorder  of  the  system”  [13],  in  other  words  it  quantifies  the  amount  of  unknown 
information  in  a  thermodynamic  system.  The  statement  of  the  second  law  varies 
between  equilibrium  thermodynamics  and  non-equilibrium  thermodynamics.  In  equi¬ 
librium  thermodynamics,  the  second  law  is  stated  in  two  parts.  The  first  part  applies 
to  reversible  processes,  while  the  second  part  may  be  applied  to  irreversible  process. 

In  equilibrium  thermodynamics,  when  a  reversible  process  affects  a  closed  sys¬ 
tem,  the  change  in  entropy  is  given  by  [39]: 


(2.29) 
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and  when  the  same  system  is  affected  by  a  irreversible  process  the  change  in  entropy 
is  given  by  [39]: 


SB-SA> 


irrev 


(2.30) 


The  first  part  explains  how  to  calculate  the  difference  in  entropy  caused  by 
a  reversible  process  changing  equilibrium  state  A  into  an  equilibrium  state  B.  The 
second  part  mandates  that  in  an  irreversible  process,  the  value  of  the  integral,  f  ^r, 
is  less  than  the  change  of  entropy  as  calculated  for  a  reversible  process. 

In  non-equilibrium  thermodynamics,  the  first  part  of  the  second  law  states  that 
the  change  in  entropy  can  be  divided  into  two  parts,  shown  in  Equation  (2.31)  [39]: 


dS  =  deS  +  diS 


(2.31) 


where  dsE  is  entropy  lost  from  the  system  to  its  surroundings  and  dtS  is  the  entropy 
produced  by  the  irreversible  process  within  the  system. 

The  second  law  also  states,  for  a  non-equilibrium  system,  the  amount  of  entropy 
produced  is  never  less  than  zero,  if  the  entropy  produced  equals  zero  the  process 
affecting  the  system  is  reversible,  and  if  the  amount  of  entropy  produced  is  greater 
than  zero  the  process  affecting  the  system  is  irreversible  [39]. 


{diS  =  0  ( reversible  process ) 
diS  >  0  {irreversible  process ) 


(2.32) 


The  discerning  reader  will  notice  the  entropy  produced  in  an  irreversible,  or  non¬ 
equilibrium,  process,  is  calculated  using  only  equilibrium  thermodynamics  in  Equation 
(2.30).  This  is  possible  because  in  equilibrium  thermodynamics,  the  entropy  produc¬ 
tion  is  found  using  the  difference  between  two  equilibrium  states.  In  non-equilibrium 
thermodynamics,  the  entropy  production  is  found  at  an  instantaneous  point  in  time. 
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The  third  part  of  the  second  law  states  that  deS  of  a  closed  system  affected  by 
either  a  reversible  or  an  irreversible  process  can  be  calculated  with  the  relationship 
[39]: 


deS 


dQ 

~T 


( closed  system,  any  process ) 


(2.33) 


The  differential  forms  of  the  second  law  in  Equations  (2.31  and  2.33)  indicate 
the  system  is  in  a  state  of  non-equilibrium,  opposed  to  the  equilibrium  form  of  the 
second  law  comparing  two  equilibrium  end  states.  Also  note,  in  non-equilibrium 
thermodynamics,  the  irreversible  statement  assumes  that  meaningful  values  of  the 
temperatures  and  entropy  exist  and  that  their  values  can  be  explicitly  solved  for.  This 
condition  limits  the  applicability  of  the  irreversible  statement  to  certain  problems  [39] . 

The  third  law  of  thermodynamics  states,  ’’The  entropy  of  any  system  in  a  state 
of  stable  equilibrium  tends  to  a  finite  value  as  the  absolute  temperature  tends  to 
zero  [34].”  The  third  law  indicates  you  don’t  need  to  know  that  absolute  values  of 
the  entropy  in  a  system  at  end  states  A  and  B,  it  is  only  the  change  in  entropy 
between  the  end  states  that  is  important.  The  third  law  also  reduces  the  amount  of 
needed  data  necessary  to  calculate  the  thermodynamic  properties  of  the  system  after 
it  has  undergone  changes.  Also,  if  the  chemical  composition  of  the  system  remains 
constant  as  the  system  undergoes  a  change,  the  third  law  allows  for  thermodynamic 
calculations  when  only  the  respective  values  of  entropy  above  a  reference  state  are 
known. 


2-4-2  Equilibrium  Thermodynamics.  According  to  equilibrium  thermody¬ 
namics  when  a  system  is  in  an  equilibrium  state,  the  state  variables  are  mathemat¬ 
ically  related  to  each  other  through  equations  of  state  (EOS).  An  EOS  consists  of 
independent  state  variables  and  dependent  state  variables,  or  state  functions.  With 
an  EOS  it  is  only  necessary  to  know  the  value  of  a  few  of  the  state  variables  in  order 
to  be  able  to  describe  the  thermodynamic  state  of  the  system  [34] . 
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A  static  system  can  be  considered  to  be  in  thermodynamic  equilibrium  when 
several  criteria  are  met.  These  criteria  consist  of  mechanical  equilibrium,  thermal 
equilibrium  and  chemical  equilibrium.  A  state  of  mechanical  equilibrium  indicates 
there  are  no  unbalanced  forces,  internally  or  externally,  in  the  system.  A  state  of 
thermal  equilibrium  indicates  the  whole  system  is  at  the  same  temperature  matching 
the  temperature  of  the  environment.  Finally,  a  state  of  chemical  equilibrium  indicates 
there  are  no  spontaneous  changes  to  the  internal  structure  of  the  system  [30].  A  dy¬ 
namic  system  may  also  be  considered  to  be  in  a  state  of  equilibrium  if  it  proceeds  slow 
enough  the  system  consists  of  multiple  stable  equilibrium  states.  After  a  shockwave 
passes,  the  value  of  the  thermodynamic  variables  will  continue  to  change  for  a  speci¬ 
fied  period  of  time  in  a  process  called  relaxation.  If  the  order  of  the  time  between  the 
equilibrium  states  is  greater  than  the  order  of  the  relaxation  time  the,  process  may 
be  considered  slow,  allowing  the  use  of  equilibrium  thermodynamics.  When  a  process 
proceeds  such  that  it  cannot  be  considered  as  a  set  of  equilibrium  states  of  a  system, 
it  is  considered  to  be  a  non-equilibrium  process.  These  equilibrium  states  generally 
do  not  happen  in  nature,  however  the  effort  of  considering  all  the  complexities  of  non- 
equilibrium  systems  may  not  be  worth  the  effort,  justifying  the  use  of  equilibrium 
thermodynamics.  [34] 

According  to  the  previous  conditions  for  equilibrium,  the  region  directly  behind 
a  shockwave,  due  to  the  relaxation  times,  cannot  be  considered  to  be  in  a  state  of 
equilibrium.  In  order  to  describe  the  thermodynamics  conditions  behind  a  shock  wave, 
however,  equilibrium  based  thermodynamics  and  some  irreversible  thermodynamic 
theories  have  introduced  the  assumption  of  a  local  equilibrium.  The  local  equilibrium 
assumption  allows  state  functions  to  be  described  in  a  ’’local”  region  using  the  same 
state  variables  necessary  to  describe  an  equilibrium  state.  Even  though  a  state  of 
local  equilibrium  is  assumed,  the  overall  system  still  produces  entropy,  the  change 
in  entropy  of  a  system  in  local-equilibrium  is  described  using  the  Gibbs  equation  for 
entropy  production: 

dS  =  ^de  +  ^ PdV  (2.34) 
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The  terms  on  the  right  hand  side  of  Equation  (2.34)  refer  to  the  contributions  of 
temperature  and  pressure  to  the  production  of  entropy  [18]. 

The  thermodynamic  state  of  a  system  can  also  be  described  using  Statistical 
Mechanics,  in  which  the  behavior  of  behavior  of  many  molecules  and  atoms  are  consid¬ 
ered.  The  collective  behavior  of  the  microscopic  particles  are  then  averaged,  producing 
statistical  concepts  of  pressure,  temperature,  internal  energy  and  entropy  [34],  For 
example  entropy  is  found  by  relating  the  distribution  of  the  microscopic  particles  to 
the  macroscopic  thermodynamics  using  Boltzmann’s  relation.  [39]: 

S  =  klnQ  (2.35) 


where  k  is  Boltzman’s  constant  and  Q  is  the  number  of  ways  to  distribute  the  particles 
in  a  system  [13].  Also,  in  statistical  mechanics  atoms  are  modeled  as  quantized 
oscillators,  with  each  oscillator  containing  three  directions  of  vibration  and  possessing 
a  mean  energy  e  [29] . 

When  a  system  is  in  either  global  equilibrium,  meaning  the  entire  system  is 
in  equilibrium,  or  local  equilibrium,  meaning  only  an  isolated  part  of  a  system  is  in 
equilibrium,  the  frequencies  atoms  oscillate  at,  geq,  are  described  by  the  Maxwell- 
Boltzman  distribution  function  [18]: 


9eq  ^ 


(2nkT\ 

\  m  J 


(2.36) 


where  n  is  the  number  of  atoms  per  unit  volume,  m  is  the  mass  of  the  atom  and  C  is 
the  relative  velocity  of  the  atom  with  respect  to  the  mean  velocity  of  the  system  [39] . 


2-4-3  Non- Equilibrium  Thermodynamics  (NET).  If  a  system  undergoes  a 
change,  and  the  rate  of  said  change  is  relatively  slow  when  compared  to  the  time 
required  for  a  system  to  relax  to  an  equilibrium  state,  the  applicability  of  the  local 
equilibrium  assumption  is  no  longer  valid  and  the  system  is  considered  to  be  in  a  state 
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of  non-equilibrium.  A  shock  wave  traveling  through  a  solid  is  an  example  of  such  a 
non-equilibrium  process  [18]. 

The  study  of  non-equilibrium  thermodynamics  (NET)  puts  a  great  emphasis 
on  the  determination  of  the  amount  of  entropy  produced  in  a  irreversible  or  non¬ 
equilibrium  system.  The  second  law  introduced  in  equation  2.31,  introduces  entropy 
as  a  system  variable  and  states  entropy  is  produced  when  a  system  undergoes  an 
irreversible  or  non-equilibrium  process.  This  law  can  also  be  written  in  terms  of 
equilibrium  and  non-equilibrium  entropy  [18]: 


(2.37) 


where  the  subscripts  i  and  f  indicate  the  initial  and  final  states  respectively,  and  Sprod 
represents  the  entropy  produced  when  the  system  goes  from  the  initial  to  the  final 
state. 

In  equilibrium  thermodynamics  entropy  is  considered  a  characteristic  state  func¬ 
tion.  This  means  it  is  able  to  describe  all  the  conceivable  thermodynamic  informa¬ 
tion  about  a  system,  if  the  entropy  is  expressed  using  the  system’s  extensive  nat¬ 
ural  variables  [19].  For  example  Callen  [8]  use  Equation  (2.34)  and  differentiate  it 
to  produce  equations  of  state  for  temperature  and  pressure,  T~1(e,V )  =  and 
T~1P(e,V)  =  |^,  respectively.  In  NET  it  is  assumed  that  non-equilibrium  entropy 
and  temperature  exist  and  that  they  have  the  potential  to  be  useful  in  the  description 
of  systems  in  non-equilibrium. 

As  discussed  in  in  Section  2.4.1  the  determination  of  the  entropy  being  produced 
differs  between  equilibrium  and  non-equilibrium  thermodynamics.  In  equilibrium 
thermodynamics,  the  production  of  entropy  in  an  irreversible  process  is  calculated 
between  two  equilibrium  states,  as  shown  in  Equation  (2.30): 
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However,  if  any  relaxation  is  occurring  and  the  order  of  the  relaxation  time  is  equal  to 
or  greater  than  the  order  of  the  time  between  equilibrium  states,  Equation  (2.30)  will 
be  unable  to  capture  the  relaxation  process.  When  this  is  the  case,  the  entropy  pro¬ 
duction  needs  to  be  described  instantaneously.  A  term  capturing  the  non-equilibrium 
relaxation  is  also  needed.  Equation  (2.31)  shown  here: 


dS  =  deS  +  diS 


accounts  for  both  conditions  necessary  to  solve  for  the  amount  of  entropy  production 
by  being  written  in  differential  form,  allowing  an  instantaneous  description  of  the  en¬ 
tropy  generation,  and  including  a  non-equilibrium  term  accounting  for  the  relaxation 
time. 

Due  to  the  local  equilibrium  assumption’s  lack  of  applicability,  both  classic  ther¬ 
modynamic  state  variables  and  non-equilibrium  state  variables  are  used  to  describe 
the  irreversible  processes  of  a  system  in  NET.  The  nature  of  these  extra  variables 
depends  on  the  system  being  studied.  Theories  of  irreversible  thermodynamics  allow 
for  the  identification  of  these  variables.  Modern  irreversible  thermodynamics  include, 
extended  irreversible  thermodynamics  (EIT)  and  internal  variable  theory  (IVT)  [19]. 

2.4.3. 1  Extended  Irreversible  Thermodynamics  (EIT).  EIT  incorpo¬ 
rates  dissipative  fluxes  with  the  use  of  the  classical  state  variables  to  describe  the 
thermodynamics  of  a  system.  Dissipative  fluxes  refer  to  irreversible  and  dissipative 
process  used  to  describe  the  interaction  of  a  material  with  its  surroundings  [26].  For 
example  Jou  et  al  [18]  [19]  use  the  heat  flux  q,  bulk  viscous  pressure  pv  and  the  devi- 
atoric  part  of  the  viscous  pressure  tensor  Pq,  with  corresponding  evolution  equations, 
to  describe  fluid  consisting  of  one  component. 

The  heat  flux  refers  to  the  flow  of  heat  into  and  out  of  a  system.  The  bulk 
viscosity  pressure  and  the  deviatoric  part  of  the  viscous  pressure  make  up  the  viscous 
pressure  tensor,  P  ,  through  the  relationship,  P  =  pvI  +  Pq.  In  this  relationship  I 
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is  identity  tensor.  The  total  pressure  tensor,  P  is  then  composed  of  the  pressure,  P, 
and  the  viscous  pressure  tensor,  P'  ,  using  the  relationship,  P  =  PI  +  P'  . 

Using  these  evolution  equations  the  linear  generalized  Gibbs  equation  for  en¬ 
tropy  production,  in  vector  notation,  becomes  [19]: 

dS  =  ©_1de  +  0_17 rdV  -  •  dq  -  ^ pvdpv  -  ^Pg  :  dPo  (2-38) 

where  ©  is  the  absolute  temperature,  7r  is  the  thermodynamic  pressure,  A,  £  and  p  are 
the  transport  coefficients  of  thermai  conductivity,  bulk  viscosity  and  shear  viscosity, 
respectively,  and  7i,  tq  and  r2  are  the  respective  relaxation  times  of  the  fluxes.  The 
right  side  of  Equation  (2.38)  is  composed  of  a  temperature  contribution,  a  pressure 
contribution,  the  contribution  of  the  heat  flux,  the  contribution  of  the  bulk  viscous 
pressure  and  the  contribution  of  the  viscous  pressure  tensor. 

The  evolution  for  the  heat  flux,  the  bulk  viscous  pressure  and  the  deviatoric 
part  of  the  stress  tensor  are  governed  using  the  following  equations  respectively  [19]: 


=  -(q  +  A  AT)  +  3"  XT2  A  •  Pg  +  3'XT2Apv 

(2.39) 

ToPv  =  ~{pv  +  £A  •  v)  +  3'(TA  •  q 

(2.40) 

r2(Po)  =  -(P[J  +  2t]v0)  +  23"vT(A0q)s 

(2.41) 

where  the  dots  above  the  variables  represent  time  derivatives  of  the  variables,  v  is  the 
barycentric  velocity  and  v0  is  the  deviatoric  traceless  part  of  the  symmetric  velocity 
gradient.  The  material  coefficients  of  3'  and  3"  are  functions  of  the  internal  energy 
and  specific  volume. 

2. 4-3. 2  Internal  Variable  Theory  (IVT).  In  internal  variable  theory, 
internal  variables  are  used  to  compliment  the  equilibrium  state  variables  in  the  de¬ 
scription  of  a  system,  for  example  Maugin  and  Muschik  [28]  define  a  constitutive 
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equation  in  the  form  of: 


cr  = 


v{x,®) 


(2.42) 


where  %  represents  the  equilibrium  state  variables  and  a  represents  any  internal 
variables  used.  Complimenting  the  classic  thermodynamic  variables  with  internal 
variables  produces  a  better  description  of  the  fast  changing  non-equilibrium  system. 
Internal  variables,  such  as  microscopic  dislocation  density,  are  used  to  describe  micro¬ 
scopic  effects  that  are  not  easily  observed.  Also,  unlike  the  macroscopic  equilibrium 
state  variables  and  EIT  fluxes,  internal  variables  cannot  be  controlled,  only  measured. 
Internal  variables  are  also  subject  to  evolution  laws,  for  example  a  in  equation  2.42 
follows  an  evolution  law  of  the  form: 


«  =  fix,  a)  +  g(x,a)x 


(2.43) 


The  use  of  IVT  indicates  the  assumption  of  a  local  equilibrium  state  [26] .  This 
allows  IVT  to  extend  beyond  equilibrium  thermodynamics,  but  requires  that  the 
thermodynamic  system  stay  near  an  equilibrium  state. 

There  are  three  main  differences  between  EIT  and  IVT.  These  differences  in¬ 
clude  the  selection  of  different  non-equilibrium  variables,  differences  in  the  evolution 
equations  of  said  variables  and  the  ability  to  control  the  non-equilibrium  state  vari¬ 
ables.  The  internal  variables  used  in  IVT  describe  microscopic  features  of  the  systems 
being  described,  while  the  flux  variables  used  in  EIT  are  invariably  independent  of 
the  system.  The  equations  describing  how  internal  variables  change  with  time,  or 
evolution  equations,  are  relaxation  equations  with  no  divergence  term  opposed  to  the 
flux  of  the  EIT  field  variables,  which  are  part  of  the  balance  equation  themselves. 
Finally,  the  flux  variables  used  in  EIT  can  be  controlled  externally,  whereas  internal 
variables  cannot  [19]. 
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2.5  Equations  of  State 

CTH  is  capable  of  using  several  different  equations  of  state  in  its  calculations, 
which  makes  this  particular  hydrocode  so  useful  in  investigating  the  differences  be¬ 
tween  different  EOS.  The  equations  of  state  being  investigated  are  the  Mie-Gruneisen 
EOS,  and  the  Sesame  EOS.  The  development  of  non-equilibrium  equations  of  state 
are  also  presented  here. 

2. 5. 0.3  Mie-Gruneisen  EOS.  The  use  of  the  Mie-Gruneisen  EOS 
begins  with  the  Hugoniot  curve.  The  Hugoniot  curve  is  an  experimentally  determined 
relationship  that  relates  two  of  the  following  variables:  pressure  (P),  speed  of  sound 
(c),  particle  velocity  (v)  and  volume  (V).  Of  the  six  potential  relationships  (P-c,  P- 
v,  P-V,  c-v,  v-V  and  c-V),  the  c-v,  P-V  and  P-w  are  particularly  useful.  In  this 
derivation  the  c-v  relationship  is  very  important.  The  graph  of  the  c-v  relationship 
typically  results  in  a  straight  line.  The  vertical  intercept  (c)  is  referred  to  as  the  bulk 
speed  of  sound  and  is  denoted  by  the  variable  Co-  The  slope  of  the  line,  S,  is  also 
an  important  quantity  to  know.  Figure  2.6  illustrates  how  Co  and  S  are  determined 
from  the  c-v  Hugoniot  [43]. 


Figure  2.6:  The  c-v  Hugoniot 
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Modern  methods  of  statically  producing  pressure  cannot  be  used  in  the  experi¬ 
mental  determination  of  the  shock  Hugoniots  for  pressures  greater  than  10  GPa.  This 
is  where  the  Mie-Gruneisen  EOS  becomes  valuable.  Meyers  [29]  presents  a  derivation 
of  the  Mie-Gruneisien  EOS,  which  is  included  here  to  give  the  reader  an  understand¬ 
ing  of  the  assumptions  made.  The  Gruneisen  parameter,  T,  is  derived  using  statistical 
mechanics.  In  statistical  mechanics,  the  total  mean  vibrational  energy  of  the  crystals 
making  up  a  material  is  equal  to  the  sum  of  the  product  of  the  number  of  energy 
levels,  n,  the  frequencies  of  the  quantized  oscillators,  u,  and  Plank’s  constant,  h,  over 
each  of  the  3N  number  of  oscillators  [29] : 


3  N 


3  N 


(2.44) 


The  summation  on  the  right  side  of  the  equation  states  the  total  mean  vibrational 
energy  is  equal  to  the  summation  of  the  mean  energies  of  the  total  oscillators,  er 
Equation(2.44)  also  indicates  the  mean  energy  of  a  harmonic  oscillator  is  equal  to 
the  product  of  the  number  of  possible  energy  levels,  the  frequencies  of  the  quantized 
oscillators  and  Plank’s  constant,  s  =  nhu  [29]. 

By  assuming  discrete  energy  levels,  the  mean  energy  of  an  oscillator  can  be 
determined.  In  statistical  mechanics,  an  energy  level  with  more  than  one  possible 
state  is  considered  to  be  degenerate.  Degeneracy  is  the  number  of  states  with  a 
common  energy  level  contained  in  a  state.  The  relative  probability  of  finding  a  system 
with  specified  levels  of  degeneracy,  g,  and  energy,  e*,  in  the  «th  level  of  the  system  is: 


(2.45) 


where  k  is  Boltzmann’s  constant.  The  absolute  probability  of  finding  a  system  with 
the  specified  energy  levels  in  its  ith  level  is  found  by  dividing  Equation  (2.45)  by  the 
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C-  j 

total  number  of  events,  JT  gje~kr ,  where  j  is  summed  over  all  energy  levels  [29]: 


(2.46) 


The  system’s  mean  energy,  s  is  found  by  summing  the  total  probability  and  the  mean 
energies  of  the  oscillators  over  i  number  of  energy  levels  [29] : 


(2.47) 


The  summation  term,  ]>T  g,(T  A" ,  is  refereed  to  as  the  partition  function  and  is 
designated  by  Q.  The  partition  function  is  the  ratio  of  the  total  number  of  molecules 
present  to  the  number  of  molecules  on  the  lowest  energy  level.  In  other  words  Q 
is  the  population  spread  among  the  available  energy  levels.  The  partition  function 
contains  the  information  necessary  to  describe  equilibrium  states,  meaning  if  Q  and 
its  temperature  dependence  is  known,  Q  may  be  used  to  calculate  the  thermodynamic 
properties  of  a  system  in  thermodynamic  equilibrium.  This  is  accomplished  by  relating 
the  statistical  description  of  the  system  (using  Q )  to  macroscopic  thermodynamics 
through  Boltzman’s  relationship,  Equation  (2.35),  to  determine  the  entropy  in  the 
system.  As  discussed  in  Section  2.4.3  entropy  is  a  characteristic  function  that  may 
be  used  to  describe  the  thermodynamic  information  of  a  system.  This  process  allows 
the  partition  function  to  act  as  a  tool  in  the  description  of  the  thermodynamics  in  a 
system  [13]. 

By  assuming  each  energy  level  contains  only  one  energy  state,  gt  becomes  one, 
meaning  the  levels  are  not  degenerate.  This  causes  the  partition  function  to  become 


[29]: 


(2.48) 


29 


Invoking  this  in  Equation(2.47),  the  value  of  the  mean  energy  becomes  [29]: 


£  = 


1  - 


hu 


his 

ekT  — 


(2.49) 


The  total  mean  energy  is  then  determined  by  summing  £  over  3N  number  of  energy 
levels  [29]: 

3  N  , 

e  =  Y.^t-  <2-5°) 

j= 1  ekT  —  1 

In  one  gram-atom  of  material,  the  total  energy  is  given  by  the  sum  of  the  vibrational 
energy,  including  the  energy  at  the  ground  state,  and  the  potential  energy  of  the 
atoms,  (j),  where  [29]: 

4>  =  -hu  (2.51) 

The  total  energy  of  one  gram-atom  of  material  becomes  [29]: 


3N  i  3N  3N 

E  =  <f>(y)  +  9 hvi  +  52  nihui  =  ^  +  52 


3= 1 


3  = 1 


L 


2/,,h  +  ~ 


3_ 

e  kT 


- 1 


(2.52) 


If  the  system  is  assumed  to  be  at  chemical  equilibrium  and  constant  tempera¬ 
ture,  the  entropy  will  tend  to  reach  a  maximum,  while  the  energy  tends  to  reach  a 
minimum.  Free  energy  is  a  state  function  that  describes  the  balance  between  these 
two  state  variables.  There  are  two  types  of  free  energy,  the  Helmholtz  free  energy 
and  the  Gibbs  Free  Energy.  If  a  constant  volume  is  assumed,  the  type  of  free  energy 
is  defined  as  the  Helmholtz  free  energy,  A,  which  is  calculated  using  the  following 
relationship  [13]: 

A  =  e  —  TS  (2.53) 

If  a  constant  pressure  is  assumed,  the  Gibbs  Free  energy  is  used  [13]: 


G  =  H  —  TS 


(2.54) 
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where  H  is  the  thermodynamic  function  of  enthalpy  equal  to  the  sum  of  internal 
energy  and  the  product  of  pressure  and  specific  volume,  H  =  e  +  PV. 

Assuming  a  constant  volume  and  inserting  the  total  energy  from  Equation  (2.52) 
into  the  Helmholtz  free  energy  produces  the  following  equation  [13]: 

E  3N  1  3N  hv 

A  =  —kT  In  =  <p(y)  +  -hzq  +  kT  In  ^1  —  e~~^  j  (2.55) 

i= 1  i=1 


By  differentiating  Equation  (2.55)  with  respect  to  volume  at  a  constant  temperature, 
the  equilibrium  value  of  pressure  can  be  found  [13]: 
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deb 

dV  + 


1  3N 
3= 1 
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hvj  + 


hu3 

hv>j 

e~kr  —  1 


(2.56) 


A  relationship  between  pressure  and  volume  has  now  been  established,  thus  complet¬ 
ing  the  relationship  between  microscopic  statistical  mechanics  and  the  thermodynam¬ 
ics  of  the  macroscopic  system.  The  term  Tj  is  [13] : 

=  _Vfdnl\  .  fdlnuA 

J  Uj\dvJT  \dlnV  J  T 

A  crystal’s  resistance  to  compression  causes  the  frequency  the  particles  are  vibrating 
at  to  increase  with  an  increase  in  the  pressure  applied  to  the  crystal.  Since  the 
vibrational  frequency  increases  as  volume  decreases,  T  >  0.  Gruneisen  assumed  all 
oscillators  in  a  system  have  the  same  T,  named  the  Gruneisen  constant  [13]: 


/  dlnu\ 

\JhPv)T 


(2.58) 


Making  all  T  equal  in  the  pressure  equation  allows  for  the  removal  of  the  summation 
term  from  Equation  (2.56),  thus  [13]: 


P  = 


dtp 

dV 


yEVIB 


(2.59) 
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When  this  pressure  equation  is  applied  atS  a  temperature  of  0  K  it  becomes  [13]: 


(2.60) 


Subtracting  Equation  (2.60)  from  Equation  (2.59)  yields  [13]: 


y(E  -  Eqk) 


(2.61) 


This  equation  is  a  form  of  the  Mie-Gruneisen  EOS  that  relates  P,  V  and  E  to  the 
pressure  and  internal  energy  at  0  K.  A  point  on  the  Hugoniot  plot  may  also  be  used 
in  which  case  the  Mie-Gruneisen  EOS  becomes  [13] : 


(2.62) 


Equations  (2.61)  and  (2.62)  can  each  be  divided  into  two  parts  [10].  The  first 
is  either  a  Hugoniot  or  a  zero-Kelvin  isotherm,  defined  in  Section  2. 5. 0.4,  which  rep¬ 
resents  the  pressure  at  a  known  state.  The  second  part  of  the  equations  is  the  con¬ 
tribution  of  vibrational,  or  thermal  energy  of  the  crystal  lattices.  In  other  words  the 
Mie-Gruneisen  EOS  can  be  represented  mathematically  as  [22] : 


(2.63) 


P  —  Pq  +  Pvib 


2. 5. 0.4  Sesame  EOS.  CTH  also  contains  the  Sesame  tabular  equation 


of  state  in  which  the  thermodynamic  state  of  the  system  is  described  using  tables 
within  CTH.  This  allows  CTH  to  reach  solutions  for  conditions  where  materials  are 
present  in  various  phases  due  to  extreme  pressures  and  temperatures.  Kerley  [20]  [21] 
provides  detailed  descriptions  on  how  the  Sesame  EOS  tables  are  built  for  iron  and 
aluminum  respectively.  The  methods  used  to  built  both  tables  are  contained  in  this 
section. 
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The  Sesame  EOS  for  iron  was  developed  by  first  constructing  the  EOS  for  the  a, 
7  and  e  phases  of  iron  using  the  PANDA  code,  which  provides  the  capability  to  model 
materials  with  complex  phase  diagrams.  These  solid  phases  include  contributions  from 
the  zero-Kelvin  isotherm,  lattice  vibrations,  and  thermal  electronic  excitations.  The 
ol  phase  also  includes  a  magnetic  contribution.  Thus  the  thermodynamic  functions  of 
pressure,  internal  energy  and  Helmholtz  free  energy  are  assumed  to  have  the  following 
forms  [20]: 

P(PiT)  =  Pc(p)  +  Pi(p,T)  +  Pm(p,T )  +  Pe(p,T),  (2.64) 

e(p,T)  =  ec(p)  +  et(p,T )  +  em(p,T)  +  ee(p,T)  +  Aeb,  (2.65) 

and 

A(p,T)  =  ec(p)  +  Ai(p,T )  +  Am(p,T)  +  Ae(p,T )  +  Ae&,  (2.66) 

Where  the  subscripts  c,  1,  m  and  e  represent  the  contributions  from  the  zero-Kelvin 
curve,  lattice  vibrations  (including  the  zero  point  term),  magnetic  excitations  (con¬ 
tribute  to  the  a  phase  only)  and  thermal  electronic  excitations,  respectively.  A Eb 
represents  cohesive  energy,  which  is  used  to  correct  for  the  zero-point  lattice  energy. 

The  zero-Kelvin  curves  are  also  referred  to  as  cold  curves,  and  account  for  ma¬ 
terials  in  ground  states  with  the  nuclei  of  the  atoms  arranged  in  perfect  crystal  lat¬ 
tices.  This  partition  of  the  EOS  accounts  for  the  cohesive  forces  leading  to  condensed 
phases  of  materials  and  the  repulsive  forces  that  determine  a  materials  response  to 
compression.  This  term  is  also  the  largest  contributor  to  the  EOS  over  much  of  the 
density-temperature  regions  commonly  studied.  [22] 

The  zero-Kelvin  curves  are  developed  by  taking  the  experimentally  determined 
compression  curves  for  the  a  and  e  phases  developed  similarly  to  the  Hugoniot  curves 
mentioned  in  Section  2. 5. 0.3,  and  subtracting  the  contributions  from  Pi  and  Pe  as 
they  are  calculated  in  the  subsequent  paragraphs.  The  remaining  pressure  is  then  fit 
to  the  Birch-Murnaghan  isothermal  EOS,  which  defines  pressure  in  terms  of  the  ratio 
of  the  density  to  the  reference  density,  rj  =  p/p0  and  the  isothermal  bulk  modulus, 
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A)  [42]: 


(2.67) 


(jl*  ~  i)  (Po  ~  4) 


For  the  7  phase,  p0  is  fixed  by  fitting  the  thermal  expansion  data,  and  the  values  of 
f3o  and  0Q  are  taken  as  intermediate  values  between  the  values  of  the  a  and  e  phases. 
Panda  is  then  used  to  compute  a  thermodynamically  consistent  energy.  [20] 


The  Lattice  vibration,  or  nuclear  degree  of  freedom  (DOF)  contribution,  is 
calculated  using  the  popular  Debye  model.  Debye  hypothesized  the  strong  bonds 
in  a  solid  dictate  the  oscillations  of  one  atom  affects  those  of  other  atoms,  and  an 
idealized  frequency  distribution  of  these  vibrations  is  required  to  model  these  effects. 
The  Debye  distribution  is  built  defining  a  maximum  frequency  vibration  of  vm  that 
acts  as  an  upper  limit,  and  then  assigns  values  to  the  rest  of  the  frequencies  [14]: 


g{y)  =  127rjK  (if  v  <  um) 
g(u)  =  0  (if  v  >  vrn) 


(2.68) 


where  Ep  is  the  potential  energy,  v  is  the  vibration  frequency,  is  an  upper  limit 
of  vibration  frequency  and  c  is  the  wave  velocity.  The  use  of  the  Debye  distribution 
allows  for  the  description  of  the  principle  features  of  the  lattice  vibration  term  by 
using  simple  select  parameters.  The  parameters  used  to  find  the  lattice  vibration 
contribution  to  the  overall  EOS  includes  the  Debye  temperature,  ©re/  =  [14], 

the  Gruneisen  parameter  at  room  temperature,  Tref,  the  density pref,  and  a  constant, 
r  which  helps  define  the  Gruneisen  function’s  dependence  on  density.  The  lattice 
contribution  to  the  Sesame  EOS  is  calculated  using  the  Sandia  code  Panda,  the  inner 
workings  of  which  are  unavailable,  however  a  simple  form  of  the  energy  contribution 
is  found  using  [22]: 


el(p,T)=Nk 


9  ©re/ 
8 


+  3 TD  (xm) 


(2.69) 
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where  D  represents  the  Debye  function  [42]: 


'm  x3dx 


(2.70) 


e*  -  1 


and  xm  =  The  pressure  contribution,  Pi  is  found  using  by  multiplying  the 

density,  Gruneisen  parameter  and  Equation  (2.69)  [22]: 


Pl{p,T)  =  U  (/>)<'! 


(2.71) 


This  assumes  the  Gruneisen  parameter  is  a  function  of  density  and  is  found  using  the 
relationship  [20]: 


(2.72) 


Note  that  when  the  EOS  accounts  for  only  the  zero-Kelvin  curve  and  the  lattice 
vibration  term,  it  can  be  shown  that  the  EOS  will  reduce  to  the  Mie-Gruneisen  EOS. 
For  example,  without  the  electronic  contribution  term  the  pressure  equation  looks 
like  [22]: 


P(p,T)  =  Pc(p)  +  P,(p,T) 


(2.73) 


which  is  similar  to  Eqn  (2.63). 

The  magnetic  contribution  to  the  Sesame  EOS,  was  handled  similarly  to  that  of 
Andrews  [3]  and  is  included  because  it  provides  a  significant  contribution  to  the  spe¬ 
cific  heat  of  a  iron  allowing  for  the  determination  of  temperature  and  internal  energy. 
The  magnetic  contribution  to  the  a  phase  of  iron  relies  heavily  on  the  noted  behavior 
of  the  Curie  temperature  of  a  iron.  The  Curie  temperature  is  the  temperature  above 
which  a  material  loses  its  ferromagnetic  properties  [33].  The  Curie  temperature  for 
iron  has  been  observed  to  be  independent  of  pressure,  which  leads  to  the  ability  to 
approximate  the  magnetic  contribution  to  the  EOS  as  independent  of  density.  This 
leads  to  the  conclusion  that  the  magnetic  contribution  to  the  pressure  term  is  negli¬ 
gible  (Pm  =  0).  Below  the  Curie  temperature,  the  heat  capacity  is  found  using  the 
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empirical  formula  [20]: 


(2.74) 


kT1-5 

CvJT)  = 

OO  -  J  ) 

where  T0  is  the  Curie  temperature  and  k  is  an  empirical  constant,  both  of  which  were 
adjusted  to  match  experimental  data.  In  a  departure  from  Andrews,  at  tempera¬ 
tures  above  the  Curie  point  Cym  =  0.  Equation  (2.74)  is  numerically  integrated  to 
determine  Em  and  Am . 

Impact  phenomena,  among  other  phenomena  (thermal  gradients,  explosions, 
etc...),  may  cause  the  material  to  become  hot  enough  the  electrons  become  very  ex¬ 
cited  and  consequently  leave  their  ground  energy  states.  This  electronic  excitation 
may  have  a  significant  effect  on  the  thermodynamic  properties  of  materials  at  high 
temperatures,  and  if  the  temperature  is  high  enough,  the  electronic  excitation  term 
may  dominate  the  EOS.  For  iron,  the  Sesame  EOS  model  predicts  that  thermal  elec¬ 
tronic  contribution  begins  to  take  effect  at  temperatures  above  500K.  Electronic  ex¬ 
citation  also  makes  a  significant  contribution  to  the  thermal  expansion,  which  agrees 
well  with  experimental  data  [20]  [22] . 

Two  separate  models  are  used  to  calculate  the  contribution  of  electronic  exci¬ 
tation  to  the  Sesame  EOS.  For  densities  above  6^0  the  INFERNO  model  of  Liber- 

cm a 

man  [25]  is  used,  and  for  densities  below  the  ionization  equilibrium  (IONEQ)  is 
used.  An  average  of  the  two  models  is  used  in  the  intermediate  densities.  [20] 

The  INFERNO  model,  models  the  potential  of  an  atom  by  taking  an  ion  sphere 
and  surrounding  it  with  a  positive  charge  and  an  electron  gas  simulating  the  neigh¬ 
boring  atoms.  Then  uses  statistical  and  quantum  mechanics  to  develop  the  desired 
EOS.  The  energies  and  wave  functions  for  the  discrete  and  continuum  electronic  levels 
are  obtained  by  solving  the  Dirac  equation  for  kinetic  energy  of  electrons  [20].  The 
INFERNO  model  simplifies  the  computation  of  the  electronic  excitation  term  at  high 
densities  [22], 

Ionization  equilibrium  theory  handles  a  system  as  a  mixture  of  free  electrons 
and  ions  at  its  simplest  form  [22].  The  free  electrons  are  treated  as  an  ideal  gas, 
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while  the  ionic  contributions  are  modeled  using  partition  functions.  The  state  of 
ionization  is  described  using  temperature  and  density,  and  is  determined  by  either 
minimizing  the  free  energy  or  balancing  the  rates  of  ionization  and  recombination 
of  the  atoms  involved.  Both  methods  assume  thermodynamic  equilibrium,  thus  the 
name  ionization  equilibrium  theory.  The  calculation  of  the  electronic  excitation  term 
using  IONEQ  is  numerically  expensive,  limiting  its  use  to  low  densities. 

A  table  of  iron’s  fluid  phase  was  built  with  the  thermal  electronic  term  and  fluid 
perturbation  theory,  which  treats  the  atomic  motions  of  the  material.  The  fluid  phase 
EOS  also  describes  the  material  in  the  vapor  and  supercritical  regimes  [20]. 

The  Sesame  tables  for  the  different  phases  of  iron  are  built  using  the  theories 
listed  above.  Then,  where  possible,  experimental  data  are  used  to  match  the  material 
parameters  used  in  iron’s  Sesame  EOS  [20]. 

A  Sesame  EOS  also  exists  for  aluminum,  in  which  the  equations  for  pressure, 
internal  energy  and  Helmholtz  free  energy  are  given  by  [21]: 

P(p,T)  =  Pc(p)  +  Pi(p,  T)  +  Pe(p,T )  (2.75) 

e(PiT)  =  ec(p)  +  ei(p,T)  +  ee(p,T)  (2.76) 

and 

A{p,T)  =  ec(p)  +  Ai(p,T )  +  Ae(p,T )  (2.77) 

The  subscript  c  stands  for  the  contributions  of  the  zero-Kelvin  curve,  represented 
by  the  analytic  formula  [21]: 

ec(p)  =  -Eo  (1  +  £  +  0.05£3)  e“*  (2.78) 

where 


37 


and  p0  is  the  density,  E0  is  the  cohesive  energy,  and  /3Q  is  the  solid  bulk  modulus,  all 
taken  at  zero  pressure. 

The  subscript  1  indicates  the  contribution  from  the  lattice  vibrations,  and  the 
subscript  e  indicates  the  contribution  of  electronic  excitation.  Both  of  which  are 
computed  in  the  same  way  the  respective  iron  contributions  are  calculated. 

The  fluid  phase  of  Aluminum  is  calculated  using  the  CRIS  [23],  model  which 
includes  contributions  from  the  ground  state  electronic  degrees  of  freedom  and  the 
contribution  of  the  nuclear  degrees  of  freedom.  The  contribution  from  the  CRIS  model 
is  added  to  the  electronic  contribution  to  describe  the  thermodynamic  state  of  the 
fluid  [21]: 

Afiuid(Pi  T)  =  An(p,  T )  +  Ae(p,  T )  (2.80) 

where  n  indicates  the  contribution  of  the  CRIS  model. 

2. 5. 0.5  Multi-Phase  Equations  of  State.  The  numerical  modeling  of 
impacts  involving  iron  present  unique  numerical  difficulties  due  to  the  fact  that  at 
pressures  of  approximately  13  MPa,  iron  undergoes  a  solid-solid  phase  change  from 
its  original  a  phase  to  the  e  phase  [5].  In  order  to  model  these  different  phases 
accurately,  two  different  sets  of  Mie-Gruneisen  parameters  are  needed.  In  CTH  you 
are  only  able  to  use  the  parameters  for  one  phase  at  a  time.  This  means  great  care 
is  needed  in  selecting  which  set  of  Mie-Gruneisen  parameters  to  use.  However,  CTH 
has  a  composite  EOS  named  Phase  Transition  Reactive  Burn  (PTRAN)  that  enables 
a  material  exhibiting  multiple  phases  to  be  modeled.  PTRAN  is  referred  to  as  a 
composite  model,  because  it  is  constructed  from  basic  EOS  models  such  as  Sesame 
and  the  Mie-Gruneisen  EOS.  The  PTRAN  EOS  uses  a  phase  boundary  to  assist  in 
describing  the  pressure  of  iron  in  the  phase  transition  region  using  the  formula  [7]: 

P(p,  T,  A  )  =  PT  +  Pt  (l  ~jj+  At  (T  ~  To)  +  AaA  (2.81) 
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where  Pt  and  px  are  the  pressure  of  the  transition  and  the  density  of  phase  1  or  a 
phase  of  iron,  A  is  the  mass  fraction  of  phase  2  or  e  phase  of  iron,  T0  denotes  room 
temperature,  (3t  in  the  bulk  modulus  in  the  transition  region  and  AT  and  A\  are  the 
derivatives  of  the  transition  pressure  with  respect  to  T  and  A.  PTRAN  has  the  option 
to  be  run  in  either  reversible  or  irreversible  mode,  though  the  documentation  does 
not  specify  the  theory  behind  either.  Unfortunately  there  is  limited  documentation 
on  this  EOS,  hindering  the  ability  to  understand  what  it  is  doing.  The  PTRAN  EOS 
is  also  capable  of  incorporating  some  non-equilibrium  thermodynamic  effects  into  an 
equation  of  state  around  a  phase  change  region  using  A  as  an  additional  internal 
variable. 

The  Sesame  equation  of  state  is  able  to  model  the  thermodynamic  properties  of 
impacts  involving  iron  much  better  as  it  has  tables  on  three  different  solid  phases  of 
iron.  The  multiphase  EOS  table  is  constructed  from  the  EOS  tables  of  the  individual 
phases.  In  a  phase  transition,  more  than  one  phase  is  present.  In  equilibrium  thermo¬ 
dynamics  when  two  phases  are  present  the  Gibbs  free  energy  of  the  two  phases  must 
be  equal.  When  only  one  phase  is  present,  it  is  the  one  with  the  lower  value  of  the 
Gibbs  free  energy  [2],  For  this  reason  in  the  iron  Sesame  table,  the  phase  exhibiting 
the  lower  value  of  Gibbs  free  energy,  at  a  given  temperature  and  pressure,  is  consid¬ 
ered  to  be  the  more  stable  phase  and  its  thermodynamic  properties  are  used  [20] .  The 
use  of  the  Gibbs  free  energy  indicates  a  constant  pressure  is  assumed  and  is  used  be¬ 
cause  the  thermodynamic  description  of  the  phase  used  is  specified  using  temperature 
and  more  importantly,  pressure.  The  Helmholtz  free  energy,  which  assumes  constant 
volume,  would  be  used  if  temperature  and  density  were  specified. 

CTH  allows  cells  in  the  mesh  to  be  divided  into  subcells,  with  each  subcell 
representing  different  materials.  This  allows  for  the  e  phase  of  the  iron  to  be  present  in 
the  same  cells  as  the  original  a  phase.  CTH  is  also  capable  of  modelling  the  different 
phases  of  iron  contained  in  a  cell  at  different  temperatures  and  pressures  [15].  To 
determine  the  cells  final  temperature  and  pressure  a  weighted  average  is  used.  This 
thermodynamic  setting  is  not  able  to  account  for  pressure  relaxation. 
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r 

fable  2.3:  EIT  Variable  Definitions 

Variable 

Definition 

q 

Heat  Flux 

Entropy 

pv 

Viscous  Pressure 

T}V 

ro 

Viscous  Pressure  Tensor 

7 r 

Thermodynamic  Pressure 

e 

Absolute  Temperature 

A 

Thermal  Conductivity  Transport  Coefficient 

c 

Bulk  Viscosity  Transport  Coefficient 

V 

Shear  Viscosity  Transport  Coefficient 

Tn 

Relaxation  Time  of  Respective  Fluxes 

2. 5. 0.6  Non- Equilibrium  Equations  of  State.  The  development  of  an 
EOS  that  can  describe  how  solid  materials  act  behind  shock  waves  is  of  great  interest 
and  the  irreversible  thermodynamic  theories  discussed  in  Section  2. 4. 3. 2  may  be  used 
to  help  model  such  phenomena.  Non-equilibrium  EOS  for  a  one  component  fluid  and 
developed  with  the  use  of  EIT  are  given  below  as  an  example  of  how  irreversible 
thermodynamic  theories  may  be  utilized  in  the  modeling  of  non-equilibrium  systems. 
For  convenience  Table  2.3  includes  the  variables  used  in  these  EOS,  first  introduced 
in  Section  2. 4. 3.1. 

In  their  review  of  EIT,  Jou  et  al  [19]  integrate  the  Gibbs  equation  in  the  form 
of  Equation  (2.38): 

dS  =  ©_1de  +  0_17r dV  —  D^-q  •  dq  —  ^-^-pvdpv  —  -^Pn  :  dP[] 

to  develop  non-equilibrium  equations  of  state  for  the  thermodynamic  pressure  and 
absolute  temperature.  The  EOS  for  O  is  found  by  integrating  Equation  (2.38)  with 
respect  to  internal  energy  [19]: 
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and  the  EOS  for  n  is  found  by  integrating  Equation  (2.38)  with  respect  to  specific 
volume  [19]: 
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The  use  of  these  expressions  assumes  that  non-equilibrium  definitions  of  O  and  7r 
exist.  Third  order  terms  in  these  equations  of  state  are  considered  to  be  negligible 
and  therefor  are  not  included.  Had  second  order  flux  terms  been  included  in  Equation 
2.38,  a  third  order  term  would  be  present. 

Another  method  of  producing  a  quasi  non-equilibrium  EOS,  is  including  an 
equilibrium  EOS,  such  as  the  Mie-Gruneisen  or  Sesame  EOS,  inside  a  rate  equation 
enabling  it  to  describe  the  thermodynamics  of  a  system  away  from  equilibrium  .  Lu 
and  Hanagud  [27]  [26]  use  this  method  in  their  works  detailed  in  the  following  section. 
Although  this  method  allows  the  EOS  to  describe  the  system  instantaneously,  the 
equilibrium  EOS  are  still  limited  by  the  equilibrium  thermodynamic  assumptions 
used  in  the  development  of  the  respective  EOS.  An  example  of  such  equilibrium 
assumptions  is  the  use  of  the  ionization  equilibrium  theory  in  the  Sesame  EOS,  which 
assumes  the  rates  of  ionization  and  recombination  occurring  are  equal  and  frozen. 
In  reality,  the  rates  of  ionization  and  recombination  are  not  equal  to  each  other  and 
change  over  time.  Thus,  this  assumption  limits  the  ability  of  the  Sesame  EOS  to 
describe  a  system  in  non-equilibrium,  even  when  used  in  a  rate  equation. 


2. 6  Modeling  of  Non- Equilibrium  Thermodynamic  Systems 

Lu  and  Hanagud  [27]  [26]  model  an  impact  loading  scenario  in  the  framework 
of  continuum  mechanics  and  non-equilibrium  thermodynamics  by  combining  the  use 
of  classical  equilibrium  state  variables  and  selected  non-equilibrium  state  variables, 
with  their  evolution  equations,  in  the  framework  of  EIT  and  IVT. 
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The  set  of  governing  equations  used  in  the  modeling  of  the  non-equilibrium 
system  requires  12  equations  consisting  of  the  conservation  equations  (Equations  2.95, 
2.96  and  2.97),  equilibrium  constitutive  equations  (Equations  2.91,  2.92  and  2.86), 
non-equilibrium  constitutive  equations  (Equations  2.89,  2.90  and  2.88),  a  kinematic 
equation  (Equation  2.84),  a  geometric  equation  (Equation  2.85)  and  some  additional 
equations  (Equations  2.93  and  2.94).  The  governing  equations  used  are  given  below 
as  an  example  of  how  non-equilibrium  thermodynamics  may  be  used  to  model  impact 
phenomena.  For  brevity,  just  the  governing  equations  previously  listed  and  some 
evolution  equations  are  included  here,  sources  [27]  [26]  may  be  consulted  to  determine 
how  supporting  variables  are  calculated. 

The  non-equilibrium  components  of  the  governing  equations  are  developed  using 
IVT  and  EIT.  IVT  is  used  to  describe  the  microscopic  dislocation  of  segments  in 
the  crystal  lattice  in  an  active  slip  system,  a.  The  orientation  of  the  slip  system 
is  described  by  the  internal  variable  tensor,  Zf-  and  is  referred  to  as  the  Schmidt 
orientation  tensor.  The  other  internal  variable  utilized  is  the  dislocation  density,  p^, 
is  a  scalar  density  describing  the  collective  length  of  dislocation  segments  in  a. 

The  rate  of  plastic  deformation,  jjPlasUc  js  a  kinematic  relationship  used  to  help 
describe  the  NET  system.  It  is  a  function  of  p%,  the  magnitude  of  the  vector  describing 
the  motion  of  the  slip  system,  or  gliding  velocity,  va,  the  magnitude  of  the  Burgers 
vector,  ba  and  It  is  determined  using  the  following  equation: 


(2.84) 


a 


To  account  for  lattice  rotations,  the  spin  vector,  uol  =  — and  the  angle 
of  Zfp  F,  =  —fiijkZfj  are  utilized.  By  combining  these  terms  the  evolution  of  Zfj  is 
determined  using  the  geometric  description: 


dtoi 


(2.85) 
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where  fl  is  the  angle  of  orientation,  and  u \  represents  the  rate  of  change  of  fl  and  is 
refered  to  as  the  spin  vector. 

As  the  material  undergoes  deformation,  the  number  of  dislocations  in  a  system 
are  going  to  change,  affecting  the  value  to  the  dislocation  density.  The  evolution  equa¬ 
tion  governing  the  change  of  ff()  is  described  using  the  following  differential  equation: 

M + hrin!  =  K«  (2.86) 

ut  uxk 

where  Ka  is  a  source  term  made  of  contribution  of  the  dislocation  nucleation  by  shock 
loading,  K%  the  dislocation  multiplication,  “  and  the  dislocation  reaction,  K? : 

Ka  =  K  +  K7  +  K7  (2-87) 


The  thermodynamic  fluxes  derived  from  EIT  include  the  dislocation  motion,  vf , 
the  non-equilibrium  stress,  cr^e  and  the  heat  flux,  qt .  The  dislocation  motion  describes 
the  motion  of  the  slip  system  and  is  made  of  the  magnitude  of  the  velocity,  va  and 
the  gliding  direction,  mf .  The  evolution  of  the  dislocation  motion  is  described  using: 


{  a  Dv^_ 
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=  0  ( otherwise ) 


(2.88) 


where  is  the  local  resolved  shear  stress,  cr“xt  is  the  external  shear  stress,  a^xt  cris 
the  critical  stress,  k%  is  the  resistance  to  dislocation  motion  and  referred  to 

as  the  substantial  derivative  of  dislocation  motion,  which  in  one-dimension  consists 
of  the  local  derivative,  ^7-  and  the  convective  derivative,  v ^ .  As  Equation  (2.88) 
indicates,  the  dislocation  motion  in  a  system  is  zero  until  a  critical  stress,  aext)Cr  is 
reached. 

The  non-equilibrium  stress,  cr™,  is  a  function  of  the  elastic  strain  rate,  and 
the  relaxation  time  of  the  viscous  process,  tv.  In  solid  metallic  materials  viscous 
processes  include  the  viscous  slip  of  the  phase  and  grain  boundaries  or  viscous  lattice 


43 


stretching.  In  liquid  metals  the  viscous  process  refers  to  the  viscous  flow  of  the 
material.  Using  the  aforemention  processes  the  evolution  of  the  non-equilibrium  stress 
is  governed  by  the  following  differential  equation: 

7b((T§e)  =  +  Vijki  Um  (2.89) 

in  which  "  represents  the  material  time  derivative. 

The  evolution  of  the  heat  flux,  qi:  is  also  found  using  a  differential  equation: 

dT 

'TqQi  =  ~Qi  ~  kijfaT.  (2-90) 

where  rq  is  the  relaxation  time  of  the  heat  flux. 

In  order  to  instantaneously  define  the  value  of  pressure,  the  equilibrium  EOS, 
pe  _  pe(p^  T),  is  found  using  the  differential  equation: 

P’  =  A(t>,T)(hj?}j+B(p,T)T  (2.91) 

This  is  an  example  of  using  an  equilibrium  EOS  in  a  rate  equation  to  give  the  NET 

hydrocode  the  ability  to  instantaneously  describe  the  thermodynamics  of  the  system. 

The  deviatoric  stress,  ofj*,  is  a  function  of  the  elastic  modulus,  Gijki,  and  the 
deviatoric  portion  of  the  deformation  tensor,  U%f.  It  is  found  using  the  rate  equation: 

*tt  =  G^l<  (2-92) 

The  total  deformation  and  stress  tensors  are  found  by  summing  their  equilibrium 
and  non-equilibrium  components,  as  shown  in  Equations  (2.93)  and  (2.94). 

°ij  =  o\j  +  (2-93) 
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and 


Uij  =  Uti  +  Unei3 


(2.94) 


Finally,  the  equilibrium  based  governing  equations  consists  of  the  familiar  con¬ 
servation  equations  [26]: 
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III.  Comparison  of  Sesame  and  Mie-Gruneisen  Based  EOS 

This  portion  of  the  numerical  investigation  compares  the  ability  of  different  EOS  to 
model  impact  phenomena.  The  EOS  being  compared  consist  of  the  Mie-Gruneisen 
EOS  for  aluminum,  the  PTRAN  EOS  for  iron,  which  includes  the  Mie-Gruneisen  EOS 
of  a  and  e  iron,  and  the  Sesame  EOS  for  both  aluminum  and  iron.  The  comparison 
of  these  EOS  will  be  performed  in  four  different  parts.  The  first  part  compares  the 
ability  of  the  Mie-Gruneisen  and  Sesame  EOS  to  model  uniaxial  impacts  of  aluminum 
over  a  range  of  impact  velocities.  The  second  part  consists  of  simulating  uniaxial  iron 
plate  impact  experiments  with  the  PTRAN  and  Sesame  EOS.  The  third  part  consists 
of  comparing  the  Sesame  EOS  to  both  the  reversible  and  irreversible  PTRAN  EOS 
to  determine  the  Sesame  EOS’s  ability  to  model  irreversible  systems.  The  fourth  and 
final  part  compares  the  answers  generated  using  all  three  EOS  to  flyer  plate  data 
gathered  by  Cinnamon  [9]  while  investigating  gouging  phenomena  for  the  HHSTT. 

Non-equilibrium  thermodynamic  theories  build  upon  equilibrium  thermody¬ 
namics.  This  requires  a  thorough  understanding  of  equilibrium  thermodynamic  theo¬ 
ries.  By  comparing  the  equilibrium  thermodynamics  contained  in  the  Mie-Gruneisen 
EOS  to  those  of  the  Sesame  EOS,  a  thorough  understanding  of  how  CTH  models  im¬ 
pact  phenomena  is  built.  Understanding  these  differences  also  assists  in  determining 
the  physics  the  currently  equilibrium  VFVC  may  fail  to  capture. 

The  uniaxial  mesh  cell  width  used  in  all  the  CTH  simulations  is  0.002  cm. 
This  size  concurs  with  what  Szmerkovsky  [35]  found  to  best  model  the  continuum 
mechanics  of  the  system.  The  uniaxial  mesh  used  in  each  simulation  is  also  designed 
to  be  wider  than  the  combined  width  of  the  flyer  and  target  plates,  with  empty  space 
next  to  the  outside  of  each  plate.  This  allows  the  outside  boundaries  of  both  plates 
to  act  as  free  surfaces  allowing  for  the  appropriate  reflection  of  the  shock  waves. 

The  equilibrium  thermodynamic  values  in  each  run  are  measured  off  the  wave 
history  generated  using  tracer  points.  In  each  impact  scenario,  the  tracer  point  used 
is  located  in  the  center  of  the  target.  This  allows  the  compression  waves  to  fully 
form  in  a  region  free  from  the  interference  and  wave  reflections  seen  at  the  edges 
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Figure  3.1:  Position  of  Captured  Thermodynamic  Values 

of  the  targets.  The  values  of  the  thermodynamic  state  variables  are  measured  at 
locations  on  the  compression  wave  well  behind  the  front  of  the  wave,  as  shown  in 
Figure  3.1.  Measuring  the  value  of  the  thermodynamic  variables  at  this  location, 
assures  the  value  is  free  from  any  numerical  noise,  or  possibly  any  non-equilibrium 
effects,  thereby  assuring  the  equilibrium  values  of  the  thermodynamic  variables  are 
obtained. 

3.1  Comparing  the  Mie-Gruneisen  and  Sesame  EOS  in  the  Uniaxial 
Impact  of  Aluminum 

3.1.1  Numerical  Setup  and  Procedure  for  Aluminum  Impacts.  The  compar¬ 
ison  of  the  Mie-Gruneisen  and  Sesame  EOS  compares  how  the  respective  EOS  model 
the  uniaxial  impact  of  aluminum.  The  material  choice  of  aluminum  is  made  for  two 
reasons.  First,  the  documentation  describing  the  development  of  these  equations  of 
state  is  publicly  available.  The  availability  of  this  documentation,  specifically  the 
Sesame  EOS,  is  necessary  to  make  sure  the  Sesame  tables  are  built  using  the  same 
theories  used  in  the  development  of  the  iron  EOS.  In  addition,  said  documentation 
is  beneficial  in  the  design  of  the  numerical  simulation  contained  in  this  investigation. 
The  second  reason  is  aluminum,  unlike  iron,  is  not  polymorphic.  This  allows  for 
a  direct  comparison  between  the  Mie-Gruneisen  EOS  and  the  Sesame  EOS  without 
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Table  3.1:  Mie-Gruneisen  Coefficients  Used  in  CTH 


Coefficient 

Po 

Go 

S 

To 

cv 

Value 

2.760  % 

5330  s 

S 

1.34 

2.16  900  ^ 

Source 

[21] 

[29] 

[29] 

[21] 

[4] 

Flyer  Target 

0.13  mm  3.327  mm 

Impact  Velocity  - ^ 

Mesh  Cell  Size  =  0.002  cm 

Figure  3.2:  CTH  Mesh  Used  to  Match  the  Experimental  Setup  of  Bekinskii  and 

Khristoforov 

the  need  to  handle  the  phase  changes  occurring  in  polymorphic  materials  under  high 
temperatures  or  pressures. 

The  form  of  the  Mie-Gruneisen  EOS  used  by  CTH  is  given  in  Equation  (3.1)  [7]: 

(  P(p,e)  =  P h(p)  +  r0p0  [e  —  eH{p)\ 

{  (3-1) 

(e(p,T)  =  eH(p)  +  Cv[T-TH(p)} 

Where  the  H  subscript  denotes  values  determined  using  the  Hugoniot  and  the  0 
subscript  denotes  the  original  values  for  T  and  p.  The  CTH  EOS  package  does  not 
contain  predefined  Mie-Gruneisen  coefficients  for  aluminum,  so  said  inputs  need  to  be 
given  to  CTH.  Table  3.1.1  contains  the  Mie-Gruneisen  coefficients  manually  specified 
inside  the  CTH  input  deck  in  SI  units.  These  coefficients  come  from  several  different 
references,  each  of  which  is  noted  in  Table  3.1.1. 

To  validate  these  coefficients,  they  are  used  in  an  impact  simulation  in  CTH 
with  an  impact  velocity  of  2020  ™ ,  matching  an  impact  experiment  performed  by 
Bekinskii  and  Khristoforov  [16].  The  mesh  used  in  this  validation  attempt  is  an 
uniaxial  recreation  of  Bekinskii  and  Khristoforov’s  flyer  plate  geometry  and  is  shown 
in  Figure  3.2.  According  to  Bekinskii  and  Khristoforov  the  pressure  behind  the  shock 
waves  should  reach  a  pressure  of  18.3  GPa.  A  tracer  point  is  placed  in  the  center 
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Impact  Velocity 


Flyer 
3.0  mm 


Target 
3.327  mm 


Mesh  Cell  Size  =  0.002  cm 

Figure  3.3:  Altered  CTH  Mesh  Used  to  Match  the  Experimental  Data  of  Bekinskii 
and  Khristoforov 

of  the  target  and  is  used  to  show  the  history  of  the  pressure  wave,  which  is  shown 
in  Figure  3.4a.  Looking  at  the  figure,  it  is  obvious  that  the  numerically  generated 
pressure  of  12.6  GPa  is  well  short  of  the  predicted  18.3  GPa.  This  inaccurate  pressure 
is  caused  by  a  mixture  of  the  artificial  viscosity  and  wave  reflection.  The  artificial 
viscosity  spreads  the  front  of  the  shock  wave  over  several  cells,  allowing  CTH  to 
handle  the  discontinuity.  However,  this  also  slows  the  development  of  the  shockwave. 
While  CTH  develops  the  shock  wave  front,  a  rarefaction  wave  from  the  flyer  passes 
through  and  stops  the  shock  wave  from  completely  forming.  To  fix  this,  the  width 
of  the  flyer  plate  is  extended  to  3.0  mm,  producing  the  final  mesh  shown  in  Figure 
3.3.  Lengthening  the  flyer  plate  increases  the  time  it  takes  for  the  rarefaction  wave 
to  reach  to  the  location  of  the  tracer  point,  allowing  CTH  to  fully  develop  the  front 
of  the  shockwave. 

The  results  generated  by  CTH  using  this  second  mesh  are  located  in  Figure  3.4b. 
The  numerical  pressure  of  18.4  GPa  compares  favorably  with  the  results  of  Bekinskii 
and  Khristoforov.  For  this  reason,  the  mesh  shown  in  Figure  3.3  is  the  mesh  used  in 
the  numerical  simulations. 

The  Johnson  Cook  constitutive  relationship  is  used  to  model  the  material  re¬ 
sponse  of  aluminum  in  the  impact  scenarios.  The  pre-defined  JC  coefficients  for 
aluminum  contained  within  CTH  are  used. 

The  comparison  of  the  Mie-Gruneisen  and  Sesame  EOS  consists  of  modeling  the 
impact  of  two  aluminum  bars  over  a  range  of  impact  velocities.  The  impact  velocities 
begin  at  an  impact  velocity  of  500  and  are  increased  in  increments  of  500  ™  until 
a  max  velocity  of  6000  —  is  reached.  The  maximum  velocity  of  6000  —  is  set  to  avoid 

J  S  J  S 
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Gauge  Pressure 


Flyer  Plate  Thickness  of  0.13  mm 


Flyer  Plate  Thickness  of  3.0  mm 


Figure  3.4:  Pressures  Generated  by  Aluminum  at  an  Impact  Velocity  of  2020  m/s 
Using  Different  Flyer  Plate  Lengths 


any  melting  that  may  occur  at  the  higher  impact  velocities  [21].  Any  melting  will 
render  the  results  generated  by  the  Mie-Gruneisen  EOS  useless,  due  to  the  EOS’s 
inability  to  handle  phase  transitions  on  its  own.  Both  EOS  will  be  used  to  model  the 
impact  at  each  impact  velocity  value.  The  results  generated  for  each  EOS  will  then 
be  compared. 


3.1.2  Results  of  Uniaxial  Impact  Simulations  with  Aluminum  Comparing  the 
Mie-Gruneisen  and  Sesame  EOS.  In  this  section,  the  results  generated  by  the 
uniaxial  impacts  of  Aluminum  are  presented  and  analyzed.  First,  the  values  of  the 
thermodynamic  variables  generated  are  plotted  against  their  respective  impact  ve¬ 
locities.  This  compares  the  ability  of  the  Mie-Gruneisen  and  Sesame  EOS  to  model 
impact  phenomena.  Then  the  thermodynamic  variables  are  plotted  against  each  other 
determining  how  the  two  EOS  handle  the  thermodynamics  of  the  system. 

Figure  3.5  shows  how  the  different  thermodynamic  variables  change  with  respect 
to  impact  velocity.  Upon  study,  Figure  3.5  shows  the  numerical  values  of  temperature 
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Density  Behind  the  Compression  Wave 


Temperature  Behind  the  Compression  Wave 


a.  Density 


Internal  Energy  Behind  the  Compression  Wave 


b.  Temperature 


Pressure  Behind  the  Compression  Wave 


c.  Internal  Energy 


d.  Pressure 


Figure  3.5:  Comparison  of  Numerical  Values  of  Density,  Temperature,  Pressure  and 
Internal  Energy  Generated  by  the  Sesame  and  Mie-Gruneisen  EOS  for  Aluminum 


and  density  generated  by  the  Mie-Gruneisen  and  Sesame  EOS  diverging  as  the  impact 
velocity  increases.  This  is  in  contrast  with  the  generated  values  of  internal  energy  and 
pressure,  which  appear  to  correlate  well  over  the  entire  range  of  impact  velocities. 
This  correlation  indicates  the  pressures  and  internal  energies  generated  for  particular 
impact  scenarios  will  be  close  in  value,  regardless  of  whether  the  Mie-Gruneisen  or 
Sesame  EOS  are  used. 

Figure  3.6  gives  the  percent  difference  between  the  thermodynamic  values  gener¬ 
ated  by  the  Sesame  and  Mie-Gruneisen  EOS  over  the  entire  range  of  impact  velocities. 
The  percent  difference  of  temperature  is  represented  by  a  dashed  line,  showing  as  the 
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impact  velocity  increases  the  percent  difference  between  the  EOS  also  increases.  This 
increase  in  difference  between  the  values  generated  with  the  different  EOS  is  reflective 
of  what  is  seen  in  Figure  3.5b,  meaning  the  generated  temperatures  correlate  well  at 
impact  velocities  below  3500  ™  with  percent  differences  of  less  than  2%,  while  at 
higher  impact  velocities  the  percent  difference  increases  at  a  faster  rate.  The  density 
percent  difference  curve,  represented  with  a  solid  line,  shows  the  difference  between 
the  values  generated  using  the  two  EOS  remains  within  remain  between  2%  and  3%, 
while  showing  a  steady  increase  over  the  entire  range  of  impact  velocities.  The  density 
curve  also  reflects  what  is  seen  in  Figure  3.5a. 

The  percent  differences  of  internal  energy  and  pressure,  represented  by  the 
dashed  dot  line  and  dotted  line  respectively,  do  not  appear  to  reflect  what  is  seen 
in  their  respective  graphs  in  Figure  3.5,  meaning  the  graphs  in  Figure  3.5  appear  to 
follow  the  same  line,  while  the  respective  percent  differences  begin  with  high  values 
and  decrease  over  the  range  of  the  impact  velocities.  These  great  percent  difference 
values  at  the  low  impact  velocities  and  the  low  percent  difference  values  at  the  high 
impact  values  is  due  to  the  range  of  internal  energy  and  pressure  values  generated. 
For  example,  the  standard  deviation,  a  measure  of  variability  representing  the  average 
difference  of  two  sets  of  data  from  the  mean  of  said  data  sets  [11],  of  the  difference 
between  the  values  of  internal  energy  generated  with  the  two  equations  of  state  is 
13.9  |%  The  standard  deviation  of  the  difference  is  4.4%  of  the  internal  energy 
(283  |j)  generated  at  an  impact  velocity  of  500  but  it  is  0.3%  of  the  internal 
energy  (4659.1  |^)  generated  at  an  impact  velocity  of  6000  ™.  The  values  of  pressure 
follow  a  similar  pattern,  though  not  as  extreme.  Bottom  line,  this  all  means  the  two 
EOS  correlate  better  at  high  impact  velocities  than  they  do  at  lower  ones. 

The  effects  the  different  EOS  have  on  the  temperature  generated  over  the  range 
of  impact  velocities  is  now  considered.  Recall  from  Section  2. 5. 0.4,  below  a  certain 
temperature,  Te,  the  Sesame  EOS  of  aluminum  is  made  up  of  two  different  parts,  the 
zero-Kelvin  curve,  pc  and  the  lattice  vibration  term,  pi.  These  two  partitions  are  the 
same  partitions  that  compose  the  Mie-Gruneisen  EOS.  This  indicates  below  Te,  the 
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Figure  3.6:  Percent  Difference  in  Values  Generated  with  the  Sesame  and  Mie- 

Gruneisen  EOS  for  Aluminum  for  Given  Impact  Velocities 

numerical  values  generated  by  each  EOS  should  correspond  with  each  other.  The 
numerical  answers  generated  by  the  two  EOS  may  even  match,  if  particular  material 
values  coefficients  are  used  by  the  Mie-Gruneisen  EOS.  Above  Te,  the  Sesame  EOS 
contains  a  component  accounting  for  the  electronic  excitation  component,  pe.  This 
suggests  at  temperatures  above  Te,  the  numerical  answers  generated  by  the  Sesame 
EOS  will  diverge  from  those  generated  by  the  Mie-Gruneisen  EOS. 

For  Aluminum,  the  value  of  Te  is  approximately  940  K  [21],  which  according  to 
Figure  3.5b,  occurs  between  the  3500  ™  and  4000  ™  impact  velocities.  Figure  3.5b 
shows,  as  predicted,  at  impact  velocities  below  4000  ™  the  answers  for  the  Mie- 
Gruneisen  and  Sesame  EOS  correspond  well  with  each  other.  It  also  shows  at  impact 
velocities  above  4000  y  the  temperatures  generated  using  the  Mie-Gruneisen  EOS 
increase  at  a  faster  rate  than  those  generated  using  the  Sesame  EOS.  This  trend  is 
also  seen  in  Figure  3.6  where  the  percent  difference  between  the  numerical  values 
of  temperature  is  less  than  2%  at  impact  velocities  below  4000  ™  while  at  impact 
velocities  equal  to  or  greater  than  4000  ™  the  value  of  percent  difference  is  greater 
than  2%.  Also  apparent  in  Figure  3.6  is  an  inflection  point  in  the  percent  difference 
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curve  located  below  the  impact  velocity  of  4000  ™  corresponding  to  the  impact  velocity 
where  the  electronic  excitation  contribution  to  the  Sesame  EOS  begins. 

This  difference  in  temperature  stems  from  the  electronic  excitation  term  in  the 
Sesame  EOS.  When  the  temperature  of  a  material  becomes  high  enough,  the  electrons 
become  excited  and  leave  their  ground  states.  During  this  process  thermal  energy  is 
absorbed,  resulting  in  the  lower  temperatures  produced  above  Te  when  using  the 
Sesame  EOS. 

Figure  3.5a  is  used  to  determine  how  the  densities  generated  with  the  different 
EOS  change  over  the  range  of  impact  velocities.  In  this  figure  the  density  generated 
by  the  Sesame  EOS  is  represented  by  a  solid  line,  the  densities  generated  by  the 
Mie-Gruneisen  EOS  are  represented  by  a  dashed  line,  and  a  third  set  of  densities 
represented  by  the  dashed/dotted  were  generated  by  translating  the  Mie-Gruneisen 
generated  densities  downward.  This  third  line  allows  an  easier  comparison  of  how 
the  Mie-Gruneisen  and  Sesame  EOS  affect  the  densities  generated  over  the  range  of 
impact  velocities. 

From  Figure  3.5a  it  is  apparent  that  the  densities  generated  begin  to  diverge 
from  each  other  right  away.  To  better  determine  the  cause  of  this,  the  densities  are 
plotted  against  the  respective  temperatures  as  shown  in  Figure  3.7.  In  this  figure  it  is 
observed  that  at  lower  temperatures,  values  of  the  densities  generated  by  the  different 
EOS  are  very  close,  then  at  a  temperature  of  just  above  420  K  the  densities  begin 
to  diverge.  Since  420  K  is  smaller  than  the  value  of  Te  for  aluminum,  it  is  indicative 
that  the  density’s  divergence  is  a  result  of  the  different  theories  used  to  generate  the 
vibrational  terms  used  in  the  Mie-Gruneisen  and  Sesame  EOS.  The  density  at  which 
this  divergence  occurs  will  be  referred  to  as  pi  and  is  about  3.2  .  This  density 

divergence  is  also  apparent  in  Figure  3.6,  showing  an  inflection  point  at  an  impact 
velocity  of  about  1500  ™,  which  produces  a  density  just  under  pi. 

The  cause  of  the  divergence  at  the  earlier  temperature  is  explained  using  Equa¬ 
tions  (2.69)  and  (2.70)  from  Section  2. 5. 0.4.  These  equations  determine  the  Debye 
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Density  Behind  the  Compression  Wave 


Figure  3.7:  Densities  Generated  Using  the  Mie-Gruneisen  and  Sesame  EOS  Plotted 
Against  Their  Corresponding  Temperatures 

function  and  lattice  vibration’s  contribution  to  the  internal  energy,  respectively  .  The 
Debye  function  is  a  parameter  used  in  Equation  (2.69).  The  key  parameter  in  the  De¬ 
bye  function  is  the  ratio  of  the  Debye  temperature  to  the  absolute  temperature, 
When  this  value  is  equal  to  one  an  inflection  point  occurs.  The  Debye  temperature 
for  aluminum  is  426  K  [14],  approximately  the  same  temperature  where  the  two  EOS 
begin  to  diverge.  This  confirms  the  divergence  is  due  to  the  different  theories  used  in 
development  of  the  lattice  vibration  contributions  and  is  not  numerical. 

In  CTH  internal  energy  and  pressure  are  state  functions  of  temperature  and 
density,  not  impact  velocity.  So,  internal  energy  and  pressure  need  to  be  plotted 
against  the  thermodynamic  state  parameters  of  temperature  and  density  before  any 
observations  are  made  about  the  thermodynamics  included  in  the  Sesame  or  Mie- 
Gruneisen  EOS. 

Figure  3.8  shows  the  values  of  pressure  and  internal  energy  generated  with 
their  respective  temperatures.  Looking  at  these  graphs  it  is  obvious  that  at  a  given 
temperature  value,  the  values  of  pressure  and  internal  energy  generated  using  the 
different  EOS  begin  to  diverge.  In  Figure  3.8  it  is  quite  apparent  the  values  of  internal 
energy  produced  by  the  two  EOS  begin  to  diverge  at  a  temperature  of  490  K,  the 
value  of  Te  for  aluminum.  For  clarification,  a  heavily  dotted  line  is  inserted  in  Figures 
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Internal  Energy  Behind  the  Compression  Wave  Close  View  of  Pressure  Behind  the  Compression  Wave 


a.  Internal  Energy  b.  Pressure 

Figure  3.8:  Comparison  of  Numerical  Values  of  Pressure  and  Internal  Energy  Gen¬ 
erated  by  the  Sesame  and  Mie-Gruneisen  EOS  Graphed  with  Respect  to  the  Corre¬ 
sponding  Temperatures 

3.8a  and  3.8b  at  the  temperature  value  of  490  K.  However,  the  temperature  where 
the  divergence  begins  to  occur  on  the  pressure  graph  appears  to  be  greater  than  the 
value  of  490  K.  To  determine  if  this  is  truly  the  case,  consider  Figure  3.9,  where  it 
is  apparent  that  up  until  Te  the  values  of  pressure  generated  with  the  different  EOS 
follow  the  same  line,  then,  as  the  temperature  reaches  Te,  the  values  of  the  pressures 
begin  to  diverge,  albeit  slowly.  This  slow  rate  of  divergence  would  indicate  that  for 
temperatures  between  940  K  and  1400  K  the  contribution  of  the  electronic  excitation 
term  to  the  Sesame  EOS  for  pressure  is  minimal  and  depending  on  the  application, 
could  be  assumed  negligible  when  modeling  impact  phenomena. 

Figure  3.10  has  the  variables  of  internal  energy  and  pressure  graphed  against 
their  respective  densities  to  illustrate  how  the  different  EOS  generate  these  values 
over  a  range  of  densities.  Upon  studying  Figure  3.10,  it  is  apparent  that  the  values 
for  both  internal  energy  and  pressure  begin  to  diverge  at  the  value  of  pi,  and  then 
begin  to  diverge  more  at  a  density  of  3.6  which  is  the  density  that  corresponds  to 
Te.  Again,  a  heavily  dotted  line  is  provided  to  easily  distinguish  the  position  of  pi  in 
Figures  3.10a  and  3.10b. 
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Close  View  of  Pressure  Behind  the  Compression  Wave 


Figure  3.9:  Magnified  View  of  Pressure  Versus  Temperature  for  the  Results  Gen¬ 
erated  with  the  Sesame  and  Mie-Gruneisen  EOS 


Internal  Energy  Behind  the  Compression  Wave 


Pressure  Behind  the  Compression  Wave 


a.  Internal  Energy 


b.  Pressure 


Figure  3.10:  Comparison  of  Numerical  Values  of  Pressure  and  Internal  Energy 

Generated  by  the  Sesame  and  Mie-Gruneisen  EOS  Graphed  with  Respect  to  the 
Corresponding  Densities 
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It  is  interesting  to  note  that  when  compared  with  respect  to  temperature  as 
in  Figure  3.8,  the  values  of  internal  energy  and  pressure  do  not  show  the  divergence 
caused  by  the  different  lattice  vibration  theories  used  in  the  two  EOS.  However,  this 
divergence  is  more  obvious  when  internal  energy  and  pressure  are  compared  using 
density  as  is  the  case  in  Figure  3.10.  This  signifies  that  both  the  Sesame  and  Mie- 
Gruneisen  EOS  are  more  sensitive  to  changes  in  pressure  than  they  are  to  changes  in 
temperature. 

3.2  Comparison  of  Sesame  and  PTRAN  EOS  Using  the  Impact  of  Iron 
Plates 

3.2.1  Procedure  and  Setup  of  Iron  Impact  Simulation.  Section  3.1  lays 
out  the  method  used  to  compare  the  Mie-Gruneisen  and  Sesame  EOS  for  a  non- 
polymorphic  material,  allowing  for  a  direct  comparison  between  the  EOS.  Iron,  as 
mentioned  in  Chapter  2,  is  a  polymorphic  material  that,  when  certain  pressures  are 
reached,  undergoes  a  solid-solid  phase  transition  from  the  a  iron  state  to  the  e  iron 
state.  Since  the  Mie-Gruneisen  EOS  cannot  handle  this  phase  change  on  its  own,  the 
composite  EOS,  PTRAN,  will  be  compared  to  the  Sesame  EOS.  Recall,  PTRAN  uses 
the  a  and  e  Mie-Gruneisen  EOS  in  their  respective  pressure  regions  and  uses  Equation 
(2.81)  (shown  here)  to  solve  for  the  pressure  in  the  region  of  the  phase  transition. 

P(p,  T,  A  )  =  PT  +  St  (l  ~jj+  At  (T  -  T0)  +  AaA 

Since  the  purpose  of  this  portion  of  the  investigation  is  to  compare  the  equilibrium 
values  of  the  thermodynamic  variables  generated  by  the  PTRAN  EOS,  as  an  equi¬ 
librium  substitute  for  the  Mie-Gruneisen  EOS,  the  reversible  option  of  the  PTRAN 
EOS  is  used. 

The  Sesame  and  PTRAN  EOS  will  simulate  experiments  involving  the  impact 
of  iron,  as  is  detailed  in  Barker  and  Hollenbach  [5].  Figure  3.11  shows  the  experi¬ 
mental  setup  used  in  their  experiments.  Figure  3.11  shows  the  flyer  plate  on  the  left, 
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Figure  3.11:  Schematic  of  Iron  Plate  Experiment  Setup  [5] 


supported  by  the  projectile  body  and  an  aluminum  support  providing  added  stability 
if  the  thickness  of  the  steel  flyer  plate  is  less  than  12  mm  thick.  The  target  plate  is 
shown  on  the  right  along  with  the  velocity  probes  measuring  the  the  impact  velocity 
protruding  from  the  target  plate,  and  the  Velocity  Interferometer  System  for  Any 
Reflector  (VISAR)  laser  system  used  to  measure  the  peak  free  surface  velocity  be¬ 
hind  the  stress  wave.  The  geometries  and  impact  velocities  of  the  experiments  being 
modeled  are  given  in  Table  3.2. 

In  the  experiment,  the  relatively  small  thicknesses  of  the  flyer  plate  and  target 
compared  to  their  diameters,  allows  the  assumption  that  any  effects  from  the  outer 
surfaces  will  not  affect  wave  propagation  in  the  centers  of  the  flyer  plate  and  target. 
This  assumption  allows  the  system  to  be  modeled  uniaxially.  Initially  the  tracer 
points  used  to  measure  the  wave  history  of  the  material  are  placed  at  the  back  of 
the  target,  matching  the  experimental  setup.  However,  interference  caused  by  shock 
waves  reflecting  off  the  back  of  the  target  stunted  the  development  of  the  shock  waves. 
To  correct  this,  the  location  of  the  tracer  point  is  moved  to  the  center  of  the  target, 
allowing  the  compression  waves  to  fully  develop. 
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Table  3.2:  Summary  of  Barker  and  Hollenbach’s  Experimental  Configurations 


Experiment 

Number 

Flyer 

Thickness 

(mm) 

Target 

Thickness 

(mm) 

Impact 

Velocity 

(?) 

1 

6.330 

6.317 

991.6 

2 

6.350 

6.312 

1150 

3 

19.14 

15.82 

997 

4 

19.14 

15.82 

1247 

5 

6.314 

6.314 

1292 

6 

6.337 

6.370 

1567 

7 

12.82 

19.14 

1292 

8 

12.81 

19.13 

1557 

9 

6.337 

6.345 

1900 

10 

6.325 

6.335 

1871 

11 

6.330 

19.07 

1887 

12 

6.327 

19.06 

1908 

13 

6.320 

6.380 

671.1 

14 

6.335 

6.375 

612.7 

19 

19.15 

15.76 

1396 

The  geometry  and  impact  velocity  of  each  experiment  is  numerically  simulated 
in  CTH  using  both  the  PTRAN  and  Sesame  EOS.  The  pre-dehned  Mie-Gruneisen 
coefficients  contained  in  CTH  are  used  for  the  modeling  of  the  a  and  e  phases  of  iron 
with  the  PTRAN  EOS.  The  Zerrili-Armstrong  constitutive  relationship  is  utilized  to 
model  iron’s  material  response  to  the  impacts  using  the  predefined  ZA  coefficients  for 
iron. 


3.2.2  Results  Generated  by  the  Numerical  Simulation  of  the  Impact  of  Iron 
Plates  . 


3. 2. 2. 1  Comparison  of  the  Equilibrium  Values.  In  this  section  CTH 
results  generated  using  the  Sesame  and  PTRAN  EOS,  in  the  numerical  simulation  of 
Barker  and  Hollenbach’s  experiments  are  presented  and  discussed.  Then,  for  further 
understanding  of  the  differences  between  the  Sesame  and  PTRAN  EOS,  additional 
impact  simulations  are  performed  and  their  results  discussed. 
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Table  3.3:  Barker  and  Hollenbach’s  Experimental  Results 


Impact 

Velocity 

(?) 

Free  Surface 
Velocity 

(?) 

Specific 

Volume 

(  —  ) 

'  kq  ' 

Peak 

Stress 

(GPa) 

612.7 

308.35 

0.1197 

12.1 

671.1 

337.58 

0.1187 

13.2 

991.6 

500.30 

0.1123 

17.3 

997.0 

505.83 

0.1122 

17.3 

1150 

581.98 

0.1108 

20.4 

1247 

633.64 

0.1100 

22.6 

1292 

655.17 

0.1095 

23.6 

1292 

657.84 

0.1096 

23.7 

1396 

711.88 

0.1086 

26.1 

1557 

791.56 

0.1072 

30.1 

1567 

792.21 

0.1071 

30.4 

1871 

949.26 

0.1047 

38.6 

1887 

NA 

0.1046 

39.1 

1900 

961.54 

0.1046 

39.6 

1908 

969.02 

0.1045 

39.8 

Barker  and  Hollenbach  used  a  VISAR  laser  to  measure  the  peak  free-surface 
velocity  behind  the  stress  wave.  Hugoniot  relationships  were  then  used  to  determine 
the  peak  stresses  and  specific  volumes  behind  the  stress  wave.  Table  3.3  [5]  contains 
their  measured  velocities,  calculated  stresses  and  calculated  specific  volumes.  In  this 
section  the  results  in  Table  3.3  are  compared  to  numerically  generated  values  using 
the  procedure  laid  out  in  Section  3.2. 

The  numerically  generated  results  using  the  Sesame  and  PTRAN  EOS  are 
graphed  with  Barker  and  Hollenbach’s  results  in  Figure  3.12.  In  Figures  3.12a  and 
3.12b  it  is  evident  that  the  numerical  answers  generated  with  the  Sesame  and  PTRAN 
EOS  both  match  the  free  surface  velocity  and  peak  stress,  which  relies  heavily  on  pres¬ 
sure,  quite  well.  However,  in  Figure  3.12b  the  values  of  specific  volume  generated  by 
the  Sesame  EOS  are  lower  than  those  generated  by  the  PTRAN  EOS  in  regions  where 
the  impact  velocity  is  high  enough  to  cause  the  e  phase  of  iron  to  be  present.  Since  the 
specific  volume  is  the  reciprocal  of  density,  these  results  indicate  the  Sesame  EOS  is 
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Velocity  Behind  the  Shockwave 


Specific  Volume  Behind  the  Shockwave 


a.  Free  Surface  Velocity  b.  Specific  Volume 


Stress  Behind  the  Shockwave 


c.  Peak  Stress 


Figure  3.12:  Comparison  of  Experimental  and  Numerical  Values  of  Free  Surface 

Velocity,  Specific  Volume  and  Peak  Stress  for  Barker  and  Hollenbach’s  Impact  Exper¬ 
iments  with  Iron 

generating  higher  values  for  density  than  the  experimental  values  and  those  generated 
by  the  e  Mie-Gruneisen  part  of  the  PTRAN  EOS. 

Table  3.4  lists  the  standard  deviation  of  the  difference  between  the  experimen¬ 
tal  results  and  those  generated  using  the  Sesame  and  PTRAN  EOS.  The  standard 
deviation  of  the  difference  between  the  Sesame  and  PTRAN  EOS  is  also  included 
in  Table  3.4.  By  comparing  the  standard  deviations  in  the  table,  the  ability  of  the 
different  equations  of  state  to  generate  answers  comparable  to  experimentally  deter¬ 
mined  values  may  be  determined.  From  the  table  it  is  apparent  the  particle  velocities 
generated  using  both  EOS  compare  equally  well  with  the  experimentally  measured 
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velocities.  However,  this  is  not  the  case  for  the  numerically  generated  values  of  specific 
volume  and  stress.  According  the  the  table,  the  values  of  specific  volume  and  stress 
generated  using  the  PTRAN  EOS  are  closer  to  the  experimentally  determined  values 
than  the  specific  volumes  and  stresses  generated  using  the  Sesame  EOS.  As  previously 
discussed  the  large  differences  in  the  numerical  specific  volumes  generated  occur  at 
impact  velocities  producing  e  iron.  This  is  due  to  different  reference  densities  used 
by  the  PTRAN  and  Sesame  EOS  for  e  iron.  The  difference  in  the  numerical  stresses 
generated  are  likely  caused  by  the  dissimilar  e  reference  densities.  Though  stress  is 
not  a  thermodynamic  variable  and  therefore  is  not  found  using  a  thermodynamic  state 
function,  the  Zerilli- Armstrong  constitutive  relationship  is  a  function  of  temperature, 
which  relies  on  the  density  for  its  generation.  Thus  the  numerical  stresses  rely  on  a 
thermodynamic  parameter  and  are  affected  indirectly  by  the  difference  in  densities. 


Table  3.4:  Comparison  of  the  Standard  Deviations  Between  the  Experimental  Val¬ 
ues,  Values  Generated  Using  the  Sesame  and  Mie-Gruneisen  EOS  for  the  Impact  of 
Iron _ 


State 

Variable 

SD  of  Sesame  EOS 
&  Experiment 

SD  of  PTRAN  EOS 
&  Experiment 

SD  of  Sesame  EOS 
&  PTRAN  EOS 

Velocity  (s) 

4.3977 

4.6126 

4.0316 

Volume 

3.0105  *10-4 

1.0207  *10-4 

2.4252  *10-4 

Stress  (GPa) 

0.2005 

0.1081 

0.2504 

The  numerical  values  of  temperature  and  internal  energy  generated  by  the  dif¬ 
ferent  EOS  can  only  be  compared  to  each  other,  because  Barker  and  Hollenbach  did 
not  provide  experimentally  determined  values  of  these  variables  for  comparison.  The 
values  of  temperature  and  internal  energy  are  plotted  against  their  respective  impact 
velocities  in  Figure  3.13.  In  Figure  3.13a  its  apparent  the  values  for  internal  energy 
generated  by  different  EOS  are  close  in  value.  However,  the  temperatures  values  gen¬ 
erated  using  different  EOS  in  Figure  3.13b  are  significantly  different.  Figure  3.13b  also 
appears  to  show  the  Sesame  EOS  generated  temperature  values  starting  to  diverge 
from  the  PTRAN  EOS  values  at  the  impact  velocity  of  1557  corresponding  to  a 
temperature  of  511.2  K.  These  results  agree  with  the  Sesame  EOS  for  iron,  which  has 
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Internal  Energy  Behind  the  Shockwave 


Temperature  Behind  the  Shockwave 


a.  Internal  Energy  b.  Temperature 

Figure  3.13:  Comparison  of  Numerical  Values  of  Temperature  and  Internal  Energy 
for  Barker  and  Hollenbach’s  Impact  Experiments  with  Iron 

a  Te  of  500  K.  However  it  is  difficult  to  determine  if  this  is  actually  the  case,  because 
the  maximum  values  of  the  temperatures  generated  by  the  PTRAN  and  Sesame  EOS 
are  638  K  and  600  K,  respectively,  which  are  not  that  much  greater  than  iron’s  Te. 

Overall  the  numerical  values  generated  using  both  the  Sesame  and  PTRAN 
EOS  compare  well  with  the  experimental  results.  The  exception  being  difference  in 
the  density  values  generated  with  the  Sesame  EOS.  The  densities  produced  by  CTH 
using  the  Sesame  EOS  in  the  aluminum  impacts,  are  different  than  those  generated 
using  the  Mie-Gruneisen  or  PTRAN  EOS.  However  in  both  cases  the  pressures  and 
internal  energies  generated  by  the  different  EOS  are  similar.  This  suggests  that  the 
’’more  accurate”  EOS  depends  on  the  application,  meaning  if  the  pressures  generated 
are  the  important  output,  either  EOS  may  be  used.  However,  if  knowing  the  density 
is  important,  the  selection  of  EOS  is  more  critical.  It  also  underlines  the  importance 
of  validating  numerically  generated  solutions  to  assure  the  correct  EOS  is  being  used 
for  the  desired  application.  For  example,  in  the  case  of  the  iron  experimental  impacts, 
the  values  of  density  generated  with  the  Mie-Gruneisen  EOS  provide  the  better  match 
to  the  experimental  values.  This  indicates  the  PTRAN  EOS  is  the  EOS  of  choice  for 
iron  impact  simulations  at  the  velocities  used  by  Barker  and  Hollenbach. 
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In  order  to  determine  more  about  the  differences  between  the  PTRAN  and 
Sesame  EOS  additional  data  points  are  needed.  To  obtain  these,  additional  runs  are 
performed  at  impact  velocities  of  400,  500,  2500,  3000,  3500,  and  4000  rj.  The  impact 
values  of  400  and  500  ™  are  impact  velocities  where  the  pressures  generated  do  not 
cause  the  a  — >  e  phase  transition.  These  values  compliment  the  experimental  impact 
velocity  of  612.7  ™  by  providing  more  data  points  where  only  the  a  phase  of  iron  is 
present.  The  geometry  used  in  the  612.7  ™  is  used  for  the  impact  velocities  of  400  ™ 
and  500  The  temperatures  produced  in  the  simulation  of  Barker  and  Hollenbach’s 
experiments  using  either  EOS  are  not  much  higher  than  iron’s  Te  value  of  500  K.  The 
impact  velocities  of  2500,  3000,  3500  and  4000  ™  are  meant  to  produce  significantly 
higher  temperatures  than  500  K,  providing  clarification  of  how  the  different  EOS  act 
at  temperatures  higher  than  Te. 

Figure  3.14  presents  the  values  of  density,  temperature,  internal  energy  and 
pressure  generated  by  the  Sesame  and  PTRAN  EOS  for  the  whole  range  of  impact 
velocities. 

There  are  several  observations  that  can  be  made  while  studying  Figure  3.14.  The 
first  is  the  similar  results  generated  by  each  EOS  for  the  values  of  pressure  and  internal 
energy.  This  is  consistent  with  what  is  seen  when  comparing  the  Mie-Gruneisen  and 
Sesame  EOS  in  aluminum  impact  tests.  In  Figure  3.14b  it  is  also  apparent  that  at 
temperatures  above  500  K  (corresponding  with  impact  velocities  greater  than  1557  ™), 
the  temperature  values  generated  by  the  Sesame  and  PTRAN  EOS  begin  to  diverge. 
As  is  the  case  in  the  Sesame  EOS  for  aluminum,  the  lower  temperature  generated  in 
the  Sesame  EOS  is  due  to  the  thermal  energy  being  absorbed  as  the  iron’s  electrons 
become  excited  and  leave  their  ground  states.  Also  apparent  in  Figure  3.14a  is  the 
difference  in  the  densities  generated  for  e  iron.  This  agrees  with  the  results  seen 
in  Figure  3.12a,  however  due  to  the  greater  range  of  density  values  in  Figure  3.14a 
it  appears  the  density  values  generated  by  the  different  EOS  correlate  better  than 
originally  thought.  It  may  be  possible  to  force  the  density  values  generated  with 
the  Mie-Gruneisen  EOS  to  match  those  generated  by  the  Sesame  EOS  for  the  e  iron 
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Internal  Energy  (KJ/kg)  Density  (kg/m3) 


a.  Density  b.  Temperature 


c.  Internal  Energy  d.  Pressure 

Figure  3.14:  Comparison  of  Numerical  Values  of  Density,  Temperature,  Pressure 

and  Internal  Energy  Generated  by  the  Sesame  and  PTRAN  EOS 
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producing  impact  velocities.  This  however,  will  cause  the  Mie-Gruneisen  EOS  to 
produce  different  pressures  and  internal  energies,  emphasizing  the  need  to  perform 
validation  runs  to  assure  the  EOS  chosen  for  use  in  numerical  simulations  is  the 
appropriate  one. 

Figures  3.14a,  3.14c  and  3.14d  also  contain  a  — >  e  phase  transition  regions  in  the 
curves  generated,  denoted  by  regions  containing  several  discontinuities.  According  to 
these  graphs,  the  impact  velocities  causing  the  phase  change  he  between  612.7  r-f, 
which  is  the  last  impact  velocity  resulting  in  the  pure  a  phase,  through  the  first 
impact  velocity  resulting  in  a  fully  developed  e  phase  at  997  ™. 

To  provide  more  information  on  how  well  the  numerical  values  of  the  thermo¬ 
dynamic  variables  correlate  with  each  other,  Figure  3.15  shows  the  percent  difference 
between  the  values  of  the  thermodynamic  variables  generated  using  the  two  EOS.  Fig¬ 
ure  3.15  shows  that  the  values  of  density,  internal  energy  and  pressure  generated  by 
the  different  EOS  correlate  well  with  each  other  throughout  the  entire  range  of  impact 
velocities,  meaning  the  percent  difference  is  never  greater  than  2%.  This  indicates 
either  the  PTRAN  or  Sesame  EOS  will  produce  acceptable  value  of  internal  energy, 
density  and  pressure  when  used  to  model  impact  phenomena.  The  percent  difference 
curve  for  temperature,  represented  by  a  dashed  line,  shows  a  great  difference  in  the 
values  of  the  temperatures  generated  at  impact  velocities  producing  e  iron.  The  per¬ 
cent  difference  then  jumps  again  after  iron’s  Te  value  of  500  K  is  reached  at  impact 
velocities  greater  than  1557  So,  if  the  values  of  temperatures  generated  in  the 
simulation  of  impacts  are  of  great  importance,  the  choice  of  EOS  is  critical.  If  high 
temperatures  are  expected  to  be  produced  in  the  impacts,  the  use  of  the  Sesame  EOS 
is  necessary,  because  it  accounts  for  the  thermal  effects  from  the  electronic  excitation 
while  the  Mie-Gruneisen  and  therefore  the  PTRAN  EOS  do  not. 

Figure  3.16b  shows  the  temperatures  generated  by  both  EOS  are  affected  by  the 
phase  transition.  The  figure  also  illustrates  the  temperatures  generated  by  the  Sesame 
EOS  in  the  e  region  are  lower  than  the  temperatures  generated  by  the  PTRAN  EOS. 
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Figure  3.15:  Percent  Difference  in  Values  Generated  with  the  Sesame  and  Mie- 

Gruneisen  EOS  of  Iron  for  Given  Impact  Velocities 

This  agrees  well  with  the  percent  difference  curve  for  temperature  in  Figure  3.15. 
There  is  no  theoretical  reason  for  the  dissimilar  values  of  temperature  generated  by 
the  two  EOS  for  the  e  phase  of  iron.  This  difference  is  most  likely  caused  by  the  use 
of  dissimilar  material  properties  values  used  by  the  different  EOS  for  the  e  iron  phase. 

As  previously  mentioned,  impact  velocity  is  not  a  thermodynamic  state  vari¬ 
able,  so  further  investigation  of  the  Sesame  and  PTRAN  EOS  requires  the  relation¬ 
ships  between  the  thermodynamic  variables  of  temperature,  density,  pressure  and 
internal  energy  be  investigated.  Figure  3.17  contains  the  values  of  density,  pressure, 
and  internal  energy  generated  by  the  Sesame  and  PTRAN  EOS  graphed  against  the 
corresponding  values  of  temperature. 

Figures  3.17a,  3.17b  and  3.17c  all  show  the  respective  values  of  density,  pressure 
and  internal  energy  for  the  different  equations  of  state  diverging  at  the  Te  value  of 
500  K,  represented  by  a  heavily  dotted  line.  The  electronic  excitation  term  included 
in  the  Sesame  EOS  accounting  for  the  difference.  In  the  density  graph  (Figure  3.17a), 
a  separate  divergence  precedes  the  divergence  at  Te ,  this  divergence  is  due  to  the 
dissimilar  densities  of  the  e  phase  of  iron  used  by  the  two  EOS.  This  differs  with 
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Figure  3.16:  Close-up  of  Temperature  Values  Generated  by  the  Sesame  and  PTRAN 
EOS  Clearly  Illustrating  the  Phase  Transition  in  Iron 

the  cause  of  the  early  divergence  of  the  densities  seen  in  aluminum,  which  occurs 
at  the  Debye  temperature,  the  material  property  at  which,  the  material  lattices  are 
vibrating  at  the  maximum  frequency  allowed  by  the  Debye  Distribution  [33].  The 
Debye  temperature  does  not  have  an  impact  on  the  iron  plate  impact  simulations 
because  neither  iron  phase  ever  cross  their  respective  Debye  temperatures.  The  Debye 
temperature  of  a  iron  in  the  Sesame  EOS  is  425  K,  but  the  maximum  temperature 
for  a  iron  seen  in  the  impact  simulations  is  343  K,  occurring  at  an  impact  velocity  of 
612.7  ™.  This  indicates  the  a  iron  phase  never  reaches  its  Debye  temperature.  The  e 
phase  of  iron  has  a  Debye  temperature  of  385  K,  which  is  reached  just  as  the  a  — >  e 
phase  transition  is  completed.  This  makes  it  difficult  to  determine  if  reaching  the 
Debye  temperature  has  any  impact  on  the  thermodynamic  properties  generated  by 
the  two  EOS.  If  it  does  have  an  impact,  it  is  difficult  to  determine  what  percentage  of 
the  difference  is  caused  by  reaching  e  iron’s  Debye  temperature  and  what  percent  of 
the  difference  is  caused  by  the  different  material  properties  used  in  the  PTRAN  and 
Sesame  EOS  for  e  iron. 


Temperature  Behind  the  Shockwave 
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Density(kg/m3) 


Density  Behind  the  Compression  Wave 


Pressure  Behind  the  Compression  Wave 


a.  Density 


b.  Pressure 


Internal  Energy  Behind  the  Compression  Wave 


c.  Internal  Energy 


Figure  3.17:  Comparison  of  Density,  Pressure  and  Internal  Energy  generated  by 

the  PTRAN  and  Sesame  EOS  with  Respect  to  Temperature 
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Internal  Energy  Behind  the  Compression  Wave 


Pressure  Behind  the  Compression  Wave 


a.  Pressure  b.  Internal  Energy 

Figure  3.18:  Comparison  of  Pressure  and  Internal  Energy  generated  by  the  PTRAN 
and  Sesame  EOS  with  Respect  to  Density 

Figure  3.18  has  the  values  of  internal  energy  and  pressure  graphed  against  the 
corresponding  values  of  density.  This  figure  shows  values  generated  for  internal  energy 
and  pressure  diverging  at  density  values  corresponding  with  the  end  of  the  a  — >  e  phase 
transition.  This  is  consistent  with  the  significant  dissimilarities  between  the  densities 
generated  with  the  different  EOS  shown  in  Figures  3.14a  and  3.17b.  Figures  3.14a 
also  show  the  density  corresponding  with  the  Debye  temperature  of  e  iron,  illustrating 
the  difficulty  in  determining  if  an  additional  divergence  is  occurring  in  this  region. 

3.3  Irreversible  Thermodynamic  Effects 

In  the  previous  sections  the  reversible  option  of  PTRAN  EOS  is  utilized  to  allow 
for  a  comparison  of  the  equilibrium  thermodynamics  included  in  the  Mie-Gruneisen 
and  Sesame  EOS.  In  this  section  the  irreversible  option  of  the  PTRAN  EOS  is  used, 
and  the  results  obtained  are  compared  to  those  generated  by  the  reversible  form  of 
the  PTRAN  and  the  Sesame  EOS.  The  purpose  of  this  investigation  is  to  use  the 
irreversible  PTRAN  EOS  to  determine  if  the  the  Sesame  EOS  is  able  to  account  for 
any  irreversibilities  due  to  phase  changes. 
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3.3.1  Procedure  for  Determining  the  Ability  of  the  Sesame  EOS  to  Generate 
Irreversibilities.  The  PTRAN  EOS  with  the  reversible  option  is  used  in  the  simula¬ 
tion  of  uniaxial  impacts  for  three  different  impact  velocities  ,  1557,  2500  and  3500 
The  impact  geometries  used  are  the  same  as  those  used  in  Section  3.2  for  the  impact 
scenarios  with  the  corresponding  impact  velocities.  Table  3.5  contains  the  respective 
geometries  for  the  different  impact  velocities.  The  only  difference  in  these  simulations 
is  the  use  of  the  irreversible  PTRAN  option,  everything  else  in  the  input  deck  remains 
the  same. 


Table  3.5:  Impact  Velocities  and  Geometries  for  Irreversible  PTRAN  Investigation 


Impact 

Flyer 

Target 

Velocity 

Thickness 

Thickness 

(?) 

(mm) 

(mm) 

1292 

6.314 

6.314 

2500 

6.327 

19.06 

2500 

6.327 

19.06 

3.3.2  Analysis  of  the  Determining  the  Ability  of  the  Sesame  EOS  to  Gener¬ 
ate  Irreversibilities.  Tables  3.6,  3.7  and  3.8  contain  the  thermodynamic  values 
of  pressure,  temperature,  internal  energy  and  density  generated  using  the  irreversible 
PTRAN  EOS,  the  reversible  PTRAN  EOS  and  the  Sesame  EOS  for  each  of  the  impact 
velocities  used.  These  values  are  measured  at  a  region  well  behind  the  front  of  the 
compression  wave  to  assure  the  respective  values  are  equilibrium  values.  Comparing 
the  thermodynamic  values  generated  behind  the  shock  determines  if  the  irreversibil¬ 
ity  effects  of  the  irreversible  PTRAN  EOS  affect  the  equilibrium  thermodynamics 
variables  generated  by  the  EOS. 

Table  3.6  contains  the  thermodynamic  values  generated  at  an  impact  velocity 
of  1292  The  values  of  internal  energy  and  pressure  generated  by  the  three  EOS 
are  all  within  1%  of  each  other,  while  the  pressures  generated  are  all  within  2%  of 
each  other.  The  temperatures  generated  with  the  different  EOS  represent  the  only 
thermodynamic  variable  where  all  the  values  generated  are  not  within  a  few  percent 


72 


or  each  other.  This  discrepancy  is  caused  by  the  temperature  generated  using  the 
irreversible  PTRAN  EOS,  which  is  5.8%  larger  than  the  temperatures  generated  using 
the  reversible  PTRAN  and  Sesame  EOS. 

From  Tables  3.7  and  3.8  it  is  ascertained  the  values  of  pressure,  internal  energy 
and  density  are  all  within  1.1%  of  each  other.  Again  the  temperatures  generated  at 
the  velocities  of  2500  and  3500  ™  are  the  only  thermodynamic  variables  with  a  span 
greater  than  2%  (8.5%  for  the  impact  velocity  of  2500  ™  and  10.1%  at  an  impact 
velocity  of  3500  ™).  In  these  cases  however,  the  rouge  temperatures  are  generated 
by  the  Sesame  EOS  and  are  due  to  the  inclusion  of  the  electronic  excitation  term 
included  in  the  Sesame  EOS. 

Unfortunately  the  cause  of  the  high  temperature  generated  by  the  irreversible 
PTRAN  EOS  cannot  be  satisfactorily  explained.  High  temperature  value  aside,  the 
irreversible  PTRAN  EOS  generates  equilibrium  thermodynamic  values  comparable  to 
the  reversible  PTRAN  EOS,  and  compares  well  with  the  Sesame  EOS  also. 


Table  3.6:  Thermodynamic  Variables  Generated  by  the  PTRAN  EOS  with  Re¬ 

versible  and  Irreversible  Options  and  Sesame  EOS  for  a  Uniaxial  Impact  with  an 
Impact  Velocity  of  1292  ™ 


EOS 

Pressure  Temperature  Internal  Density 

Energy 

(GPa)  (K)  (M)  ($) 

PTRAN  irreversible 
PTRAN  reversible 
Sesame 

23.75  451.6  361.96  9.121 

23.82  425.2  361.25  9.123 

23.43  425  364.5  9.183 

One  of  the  purposes  of  this  thesis  is  to  determine  if  the  Sesame  EOS  is  capable 
of  modeling  any  irreversible  thermodynamics  whatsoever.  According  to  the  results 
generated  in  the  simulations  of  uniaxial  impacts  involving  aluminum,  it  is  apparent  no 
irreversibilities  are  included  in  the  Sesame  EOS.  However,  by  modeling  a  polymorphic 
material,  such  as  iron,  an  irreversibility  is  introduced  in  the  phase  transition.  Looking 
at  the  equilibrium  values  generated  by  the  reversible  and  irreversible  PTRAN  EOS 
and  the  Sesame  EOS,  no  conclusions  about  the  irreversible  nature  of  the  Sesame  EOS 
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Table  3.7:  Thermodynamic  Variables  Generated  by  the  PTRAN  EOS  with  Re¬ 

versible  and  Irreversible  Options  and  Sesame  EOS  for  a  Uniaxial  Impact  with  an 
Impact  Velocity  of  2500  ™ 


EOS 

Pressure  Temperature  Internal  Density 

Energy 

(GPa)  (K)  (M)  (5) 

PTRAN  irreversible 
PTRAN  reversible 
Sesame 

57.69  953.8  900  9.976 

57.5  952.4  891.67  9.963 

57.22  897.4  897.12  10.023 

Table  3.8:  Thermodynamic  Variables  Generated  by  the  PTRAN  EOS  with  Re¬ 

versible  and  Irreversible  Options  and  Sesame  EOS  for  a  Uniaxial  Impact  with  an 
Impact  Velocity  of  2500  ™ 


EOS 

Pressure  Temperature  Internal  Density 

Energy 

(GPa)  (K)  (M)  (^) 

PTRAN  irreversible 
PTRAN  reversible 
Sesame 

93  1750  1625  10.611 

92.35  1750.4  1619  10.610 

92  1590  1625  10.656 

may  be  made.  So,  now  the  ability  of  the  different  EOS  to  model  the  history  of  a 
compression  wave  is  investigated.  Figure  3.19  shows  the  wave  histories  generated 
using  the  three  EOS  discussed  in  this  section. 

Figures  3.19a  and  3.19b  show  the  irreversible  and  reversible  forms  of  the  PTRAN 
EOS.  Comparing  these  two  figures  it  is  noticed  the  front  portion  of  the  compression 
wave  is  not  affected  by  the  irreversibility  option.  The  back  of  the  compression  waves 
are  formed  very  different,  however.  The  back  of  the  compression  wave  generated 
with  the  irreversible  option  on  shows  a  rounded  portion  where  the  e  — ►  a  phase 
change  takes  place  followed  by  a  sloped  portion  formed  by  the  rarefaction  wave.  This 
contrasts  with  the  compression  wave  generated  using  the  reversible  PTRAN  EOS 
showing  a  flat  e  — »  a  region  and  a  sharp  drop  off  in  the  region  formed  after  the  phase 
change  formed  by  a  rarefaction  wave. 


74 


PRESSURE  (GPd) 


a.  Irreversible  PTRAN  EOS  b.  Reversible  PTRAN  EOS 


c.  Sesame  EOS 

Figure  3.19:  Comparison  of  Sesame  EOS  to  Reversible  and  Irreversible  Forms  of 

the  PTRAN  EOS 
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Flyer  Plate 


Containment 

chamber 

Vacuum 


Manganin  Gauge  Velocity  Pins 


Figure  3.20:  Schematic  of  Flyer  Plate  Test  Setup 


The  compression  wave  history  generated  with  the  Sesame  EOS  has  the  same 
drop  off  in  the  region  behind  the  e  — *  a  phase  transition  seen  on  the  compression 
wave  generated  using  the  reversible  PTRAN  EOS.  However,  in  the  region  of  the 
e  — ►  a  transition  it  is  rounded,  similar  to  that  of  the  irreversible  PTRAN  EOS  gen¬ 
erated  compression  wave.  This  signifies  that  the  Sesame  EOS,  in  spite  of  using  an 
equilibrium  approach  in  determining  the  thermodynamics  of  the  phase  transition  re¬ 
gion  by  selecting  the  phase  with  the  lower  Gibbs  free  energy,  is  capable  of  modeling 
irreversible  thermodynamics  in  the  region  of  a  phase  change. 


3-4  Flyer  Plate  Experiment 

3-4-1  Procedure  and  Experimental  Setup.  While  Cinnamon  [9]  was  conduct¬ 
ing  his  research  on  gouging  for  the  HHSTT,  flyer  plate  tests  were  performed  at  the 
University  of  Dayton  Research  Institute  (UDRI)  to  assist  in  the  determination  of  the 
Zerrili- Armstrong  (ZA)  coefficients  used  to  model  the  material  response  of  1080  steel. 
A  schematic  of  the  test  setup  is  shown  in  Figure  3.20,  which  shows  the  flyer  plate 
on  the  right  attached  to  a  sabot  for  support.  The  target  plate  is  on  the  left  with  a 
stress  gauge  between  it  and  a  layer  of  PMMA.  The  impact  velocity  is  measured  using 
velocity  pins  located  just  to  the  right  of  the  target  and  are  shown  protruding  through 
the  barrel  of  the  gas  gun.  To  determine  the  correct  ZA  coefficients  for  1080  steel, 
Cinnamon  produced  a  uniaxial  flyer  plate  model  in  CTH  using  the  Sesame  EOS  of 
iron  and  manually  specifying  the  ZA  coefficients  for  1080  steel.  He  then  ran  the  simu- 
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lations  while  modifying  the  values  of  the  ZA  coefficients,  until  the  numerical  answers 
produced  by  CTH  matched  the  experimental  data. 

In  this  investigation,  the  CTH  model  used  by  Cinnamon  is  run  using  the  a  iron 
Mie-Gruneisen  EOS,  the  reversible  and  irreversible  forms  of  the  PTRAN  EOS  and 
the  Sesame  EOS.  For  alloy  metals  such  as  steel,  it  is  an  accepted  approximation  to 
model  the  alloy  using  the  EOS  for  it’s  prime  constituent,  in  this  case  iron  [43].  This 
corresponds  to  what  Cinnamon  did  in  his  simulation,  and  is  the  method  used  in  this 
investigation.  Due  to  the  relatively  low  impact  velocity,  only  the  a  iron  Mie-Gruneisen 
EOS  coefficients  are  used  to  simulate  1080  steel  impacts. 

The  ZA  coefficients  for  1080  steel  used  in  this  simulation  are  located  in  Table 
3.9.  These  are  the  coefficients  determined  by  Cinnamon  to  accurately  model  the 
reaction  of  1080  steel  when  subjected  to  high  velocity  impacts. 


Table  3.9:  Zerilli-Armstrong  Coefficients  for  1080  Steel 


Coefficient 

A 

Cl 

c2 

c3 

c4 

^5 

n 

Value 

0.825 

4.0 

0 

160.0 

12.0 

0.266 

0.089 

The  different  EOS  will  be  used  to  simulate  the  impact  of  1080  steel  at  an 
impact  velocity  of  891  ™ .  The  results  obtained  from  these  simulation  will  then  be 
compared.  The  tracer  point  used  to  generate  the  stress  wave  histories  is  located  in 
the  target  6.25  mm  away  from  the  point  of  impact.  To  avoid  any  interference  from 
wave  reflections  at  the  iron/PMMA  interface,  the  PMMA  will  be  modeled  as  iron. 
This  also  matches  the  methodology  used  by  Cinnamon,  and  is  justified  because  the 
PMMA  is  located  behind  the  tracer  point  and  has  no  effect  on  the  development  of 
the  stress  wave. 

3-4-2  Results  of  Flyer  Plate  Experiment.  Figure  3.21  contains  the  stress 
waves  generated  in  the  CTH  simulation  of  the  flyer  plate  experiments  performed  by 
Cinnamon,  using  the  a  Mie-Gruneisen,  reversible  and  irreversible  PTRAN  EOS  and 
the  Sesame  EOS.  The  stress  wave  measured  in  the  lab  is  located  in  Figure  3.22, 
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superimposed  on  the  CTH  generation  of  the  stress  wave  Cinnamon  used  to  determine 
the  correct  ZA  coefficients  for  1080  Steel  [9]. 

Upon  study,  all  the  EOS  correctly  generate  the  value  of  the  Hugoniot  Elastic 
Limit  (HEL)  and  are  able  to  model  the  elastic  precursor.  This,  however  is  the  only 
similarity  between  the  numerical  results  generated  with  the  Sesame,  PTRAN  and  a 
Mie-Gruneisen  EOS.  There  are  two  reasons  the  single  phase  Mie-Gruneisen  EOS  is 
able  to  model  accurately  model  the  value  of  the  HEL.  The  first  reason  is  the  HEL 
is  a  material  response  modeled  by  the  constitutive  equations  and  therefore  is  not 
determined  by  the  EOS.  The  other  reason  is  the  value  of  the  HEL  is  well  within  the 
limits  of  the  a  phase  of  iron. 

The  differences  in  the  results  generated  using  the  Mie-Gruneisen  EOS  as  opposed 
to  those  generated  by  the  Sesame  and  reversible  PTRAN  EOS  will  be  examined  first. 
The  two  major  differences  in  the  Mie-Gruneisen  generated  results  as  opposed  to  those 
generated  with  the  Sesame  or  PTRAN  EOS  include,  the  higher  value  of  peak  stress 
generated  with  the  Mie-Gruneisen  EOS  and  the  relatively  small  depth  to  the  stress 
wave  generated. 

The  lack  of  depth  in  the  stress  wave  is  caused  by  a  rarefaction  wave  relieving 
the  stress  wave  sooner  with  the  Mie-Gruneison  generated  results  than  the  Sesame 
or  PTRAN  results.  The  velocity  of  a  shock  wave  propagating  through  a  material  is 
found  using  the  linear  equation,  Eqn  3.2  [38]: 

US  =  C0  +  Sv  (3.2) 

where  Us  is  the  shock  wave  velocity.  Us  is  linearly  related  to  the  velocity  of  the 
material  particles  located  behind  the  shock  wave.  Table  3.10  gives  the  values  of  the 
peak  particle  velocity  found  with  each  EOS.  Table  3.10  shows  the  peak  particle  speed 
generated  with  the  Mie-Gruneisen  is  approximately  32%  higher  than  those  generated 
with  the  other  EOS.  This  leads  to  a  higher  shock  wave  velocity  causing  the  rarefaction 
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a.  Reversible  PTRAN  EOS 


b.  Irreversible  PTRAN  EOS 


T1MF  (/as) 


T1MF  ()as) 


c.  Sesame  EOS 


d.  Mie-Gruneisen  EOS  for  the  a  Iron  Phase 


Figure  3.21:  Comparison  of  Stress  Waves  Generated  with  the  Mie-Gruneisen, 

PTRAN  and  Sesame  EOS 
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14.0 


Figure  3.22:  UDRI  Results  Superimposed  on  Cinnamon’s  CTH  Generated  Stress 
Wave 

wave  to  reach  the  location  of  the  tracer  point  in  less  time  resulting  in  the  smaller  depth 
of  the  Mie-Gruneisen  generated  stress  wave. 


Table  3.10: 


!3eak  Particle  Velocities  Generated  by  Different  EOS 


EOS 

Velocity  (*) 

Mie-Gruneisen 

738.5 

Sesame 

334.4 

Reversible  PTRAN 

331.7 

Irreversible  PTRAN 

329.4 

The  difference  in  peak  stress  value  is  caused  by  the  inability  of  the  Mie-Gruneisen 
EOS  to  handle  the  phase  change  of  iron.  Equation  (2.15)  in  Section  2.2.2  shows  that 
the  hydrostatic  pressure  produced  in  an  impact  is  a  major  contributor  to  the  overall 
value  of  stress  seen  by  the  material.  The  peak  values  of  the  pressures  generated  with 
each  EOS  are  given  in  Table  3.11.  The  values  given  in  Table  3.11  show  that  the 
pressures  generated  with  the  Sesame  and  PTRAN  EOS  are  approximately  equal  to 
the  transition  pressure  of  iron.  The  Mie-Gruneisen  EOS  does  not  account  for  iron’s 


80 


phase  transition  and  generates  a  higher  value  of  peak  pressure,  resulting  in  a  higher 
peak  stress  value. 


EOS 

Pressure  (GPa) 

Mie-Gruneisen 

17.6 

Sesame 

12.9 

Reversible  PTRAN 

12.8 

Reversible  PTRAN 

12.7 

Table  3.11:  Peak  Pressures  Generated  by  Different  EOS 

Comparing  the  stress  waves  generated  using  the  reversible  and  irreversible  forms 
of  the  PTRAN  EOS  shown  in  Figures  3.21a  and  3.21b  respectively  it  is  apparent  the 
reversible  -  irreversible  option  in  the  PTRAN  EOS  may  have  a  major  impact  on  the 
results  generated.  The  reversible  form  of  the  PTRAN  EOS  generates  a  stress  wave 
with  more  depth  than  the  stress  wave  produced  using  the  irreversible  form.  Unlike 
the  difference  seen  between  the  Mie-Gruneisen  and  Sesame  and  reversible  PTRAN 
EOS  generated  stress  waves,  there  is  no  obvious  physical  reason  for  the  different 
shape  of  the  stress  wave  generated  by  the  irreversible  PTRAN  EOS.  Meaning  the 
velocities  behind  the  stress  waves  generated  using  the  different  forms  of  the  PTRAN 
EOS  compare  very  well  (331.7  ™  generated  with  the  reversible  form  compared  to 
329.4  ™  produced  by  the  irreversible  form),  so  the  plastic  wave  velocities  should  also 
compare  well  by  generating  stress  waves  of  equal  depth.  In  spite  of  the  different  stress 
wave  shapes,  the  magnitude  of  the  peak  stress  generated  with  the  irreversible  PTRAN 
EOS  is  comparable  to  that  generated  using  the  reversible  form  of  the  PTRAN  EOS, 
12.7  and  12.8  GPa  respectively.  The  results  generated  using  the  two  PTRAN  EOS 
indicate  care  must  be  used  in  the  use  of  the  PTRAN  EOS  in  the  modeling  of  impact 
events  causing  phase  transitions. 

Finally,  it  has  been  determined  the  Sesame  and  reversible  PTRAN  EOS  appear 
to  generate  similiar  results.  To  verify  this,  the  results  generated  using  the  two  EOS 
are  compared  to  the  experimental  results  in  Figure  3.22.  Looking  at  the  different 
results,  it  is  apparent  that  both  the  Sesame  and  reversible  PTRAN  EOS  generate 
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reasonable  stress  wave  histories,  both  in  peak  stress  value,  just  below  13  GPA,  and 
correct  wave  depth,  around  2  jis. 
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IV.  Validation  of  the  Vanderhyde  Finite  Volume  Code 

(VFVC) 

4-1  Validation  of  VFVC  Using  Aluminum 

4-1.1  Procedure.  The  validation  of  the  VFVC  consists  of  simulating  the 
impact  of  aluminum  using  the  VFVC  and  comparing  the  answers  generated  with 
answers  generated  by  CTH  using  both  the  Mie-Gruneisen  and  Sesame  EOS.  The 
simulation  consist  simulating  the  same  set  of  impact  velocities  and  geometries  used 
in  Section  3.1.  A  cell  width  of  0.0018  cm  is  used  in  the  validation  using  aluminum. 
The  VFVC  also  requires  additional  material  properties  in  order  to  model  aluminum’s 
constitutive  response  to  impact.  These  additional  parameters  are  given  in  Table  4.1, 
where  Y  is  the  value  of  aluminum’s  yield  stress,  eo  is  the  initial  value  of  the  specific 
internal  energy  of  iron,  G  is  the  shear  modulus  of  aluminum  and  H  is  the  linear  work 
hardening  rate. 

Table  4.1:  Additional  Material  Parameters  Needed  by  the  Eulerian  Code 


Coefficient 

Y 

e0 

G 

H 

Value 

34.747  MPa 

268125  f 

kg 

25  GPa 

120  MPa 

Source 

[12] 

[21] 

[4] 

[31] 

For  every  simulation,  the  VFVC  produces  two  output  files.  One  shows  what  is 
happening  over  the  region  of  the  entire  grid,  while  the  other  output  file  provides  the 
history  of  a  single  point  over  the  run  time  of  the  simulation.  A  post  processing  tool 
written  in  Matlab  then  reads  the  output  files  and  presents  the  data  in  an  intelligible 
manner. 

The  results  generated  by  the  VFVC  are  then  compared  to  the  numerical  data 
generated  by  CTH  using  both  the  Sesame  and  Mie-Gruneisen  EOS. 

4-1.2  Validation  of  the  Finite  Volume  Code  with  Impacts  Between  Aluminum 
Plates.  Figure  4.1  shows  the  thermodynamic  properties  behind  the  shockwave 
generated  by  the  VFVC  compared  to  those  generated  with  the  Mie-Gruneisen  and 
Sesame  EOS  in  CTH.  Figure  4.1a  contains  an  additional  curve  that  consists  of  the 
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Internal  Energy  Behind  the  Compression  Wave 


a.  Density  b.  Internal  Energy 


Pressure  Behind  the  Compression  Wave 


c.  Pressure 


Figure  4.1:  Comparison  of  Density,  Internal  Energy  and  Pressure  generated  by  the 
Finite  Volume  Code  and  CTH 

density  values  generated  by  the  VFVC  translated  to  the  right  for  an  easier  comparison 
the  densities  generated  with  the  VFVC  and  CTH,  using  the  Sesame  EOS,  over  the 
range  of  impact  velocities.  Upon  examination  of  Figures  4.1a  and  4.1c,  it  is  apparent 
that  the  values  of  density,  pressure  and  internal  energy  generated  by  the  VFVC  match 
the  values  produced  by  CTH  quite  well. 

For  a  quantitative  comparison  of  how  the  well  the  VFVC  matches  the  ther¬ 
modynamic  values  generated  with  the  different  EOS  in  CTH,  Table  4.2  contains  the 
standard  deviations  (SD)  between  the  the  results  generated  with  the  VFVC,  Sesame 
EOS  and  Mie-Gruneisen  EOS  in  CTH.  The  SD  is  a  measure  of  variability  representing 
the  average  difference  of  a  sets  of  data  from  the  mean  of  said  data  sets  [11]  making 
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it  an  ideal  value  in  determining  how  the  values  generated  by  the  VFVC  and  CTH 
compare.  The  SD  in  Table  4.2  of  the  densities  and  pressures  generated,  clearly  shows 
that  the  VFVC  and  Mie-Gruneisen  EOS  match  very  well  with  the  respective  SD  of 
0.0068  kg/kJ  and  0.1471  GPa.  The  SD  of  40.257  kJ/kg  for  internal  energy  is  also 
fairly  small  considering  the  the  orders  of  magnitude  for  internal  energy  are  103  and 
104.  Table  4.2  also  shows  that  the  SD  between  the  answers  generated  with  the  VFVC 
and  CTH  using  the  Sesame  EOS  are  comparable  to  the  ST  between  the  answers 
generated  using  the  Sesame  and  Mie-Gruneisen  EOS  in  CTH.  The  values  of  the  SD 
between  the  different  data  sets  generated  lead  to  the  conclusion  that  when  modeling 
a  non-polymorphic  material  the  thermodynamic  variables  generated  by  the  VFVC 
compare  quite  well  to  the  values  generated  by  CTH. 


Table  4.2:  Comparison  of  the  Standard  Deviations  Between  the  VFVC,  Sesame 

EOS  and  Mie-Gruneisen  EOS  Generated  Values _ 


State 

Variable 

SD  of  VFVC  &  SD  of  VFVC  &  STD  of  Sesame  & 

Sesame  EOS  Mie-Gruneisen  EOS  Mie-Gruneisen  EOS 

Pressure  (GPa) 
Density  (^) 

Internal 

Energy 

0.3347  0.1471  0.3933 

0.0102  0.0068  0.0163 

39.191  40.257  44.303 

Figure  4.2  shows  the  different  wave  profiles  generated  for  pressure  at  an  impact 
velocity  of  500  ™  by  CTH  using  both  the  Mie-Gruneisen  EOS  and  the  Sesame  EOS, 
and  the  results  generated  by  the  VFVC.  Upon  inspection  the  firgure  shows  that  the 
results  generated  by  the  VFVC  behind  the  plastic  wave  are  comparable  to  those 
generated  by  CTH.  The  graphs  also  show  that  the  VFVC  fails  to  generate  an  elastic 
precursor.  This  failure  is  caused  by  the  1st  spatial  order  Lax-Friedrichs  flux  scheme, 
which  overloads  the  discontinuity  with  dissipation  causing  the  compression  wave  to 
spread  itself  out  over  several  cells  in  the  x-direction,  thereby  smearing  out  the  elastic 
precursor.  Even  with  the  addition  of  the  MUSCL  scheme,  which  produces  a  2nd  order 
spatially  accurate  result  is  not  enough  to  produce  the  elastic  precursor.  To  remedy 
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this  problem  either  a  less  dissipative  1st  order  flux  scheme  needs  to  be  used,  or  a  third 
order  scheme  needs  to  be  introduced  into  the  VFVC. 

Overall  several  main  observations  are  made  about  the  similiarities  and  differ¬ 
ences  in  the  Mie-Gruneisen,  PTRAN  and  Sesame  EOS.  This  includes  the  similar 
pressures  and  internal  energies  generated  by  the  different  EOS,  allowing  them  to  be 
used  in  place  of  each  other  if  the  value  of  pressure  or  internal  energy  is  of  primary 
interest.  This  also  signifies  the  form  of  the  Mie-Gruneisen  used  in  the  VFVC  will 
produce  usable  results.  Also,  both  EOS  are  more  sensitive  to  changes  in  density  than 
to  changes  in  pressure.  It  is  also  shown  that  the  densities  generated  with  each  EOS 
are  different  in  value,  this  may  be  corrected  by  using  different  Mie-Gruneisen  EOS 
coefficients,  this  however,  will  affect  the  pressures  produced  by  the  Mie-Gruneisen 
EOS. 

4-2  Validation  of  VFVC  for  Iron  Impacts 

4-2.1  Procedure.  The  experimental  geometries  and  impact  velocities  used 
in  Section  3.2  are  used  in  the  validation  of  the  Eulerian  Code.  In  its  current  state 
the  Eulerian  code  is  only  capable  of  handling  one  phase  at  a  time,  meaning  it  cannot 
accurately  handle  the  phase  change  of  iron.  However  the  experimental  data  being  used 
to  validate  the  code  only  has  two  impact  scenarios  that  cause  the  pressures  generated 
to  fall  in  the  phase  change  region,  and  both  of  these  impact  velocities,  671.1  and  991.6 
y,  are  located  near  the  either  beginning  or  the  end  of  the  phase  change  region.  This 
proximity  to  these  phase  change  boundaries  allows  the  VFVC  code  to  model  these 
velocities  scenarios  as  single  phases,  in  this  case,  the  a  and  e  phases  respectively. 

The  material  coefficients  of  both  the  a  and  e  phases  used  by  the  Mie-Gruneisen 
EOS  are  taken  from  Boucheron  et  al.  and  are  located  in  Table  4.3.  The  VFVC  utilizes 
a  Fortran  ”if”  statement  to  determine  which  coefficients  to  use.  This  determination 
uses  the  magnitude  of  the  pressure  generated  by  the  impact.  In  order  to  allow  the 
boundary  impact  scenarios  to  be  modeled  as  single  phases,  the  VFVC  uses  a  pressure 
value  to  15  GPa  to  differentiate  between  the  two  phases  of  iron.  This  is  done  instead 
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Figure  4.2:  Comparison  of  Pressure  waves  generated  by  the  Finite  Volume  Code 

and  CTH 
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Table  4.3: 


Mie-Gruneisen  Coefficients  of  the  a  and  e  Phases  of  Iron 


Coefficient 

“or 

Po 

S 

To 


4600  f 
7870  \ 
1.46 
1.7 


4600  f 
8290  \ 

ma 

1.51 

2.4 


of  beginning  the  phase  transition  at  a  pressure  of  13  GPa  and  completing  it  around 
a  pressure  of  17  GPa,  as  would  be  done  in  a  code  that  can  handle  two  phases. 

Material  coefficients  are  also  required  to  describe  iron’s  constitutive  response  to 
the  impact  phenomena.  These  coefficients  are  given  in  Table  4.4. 


Table  4.4:  Additional  Materials  Parameters  Needed  by  the  Eulerian  Code 


Coefficient 

Y 

e0 

G 

H 

Value 

50.0  MPa 

131800  f 

kg 

87.0  GPa 

2600  MPa 

Source 

[4] 

[21] 

[4] 

Extrapolated 
from  [32] 

After  all  the  experimental  geometries  and  velocities  are  completed,  the  results 
generated  by  the  VFVC  are  compared  to  the  results  numerically  generated  by  CTH 
using  both  the  PTRAN  and  Sesame  EOS. 

4-2.2  Validation  of  the  Finite  Volume  Code  with  Impacts  Between  Iron  Plates. 

Figure  4.3  shows  the  thermodynamic  variables  of  density,  internal  energy  and  pres¬ 
sure  plotted  against  the  impact  velocities  producing  said  thermodynamic  variables. 
The  densities  generated  by  the  different  codes  and  EOS  are  located  in  Figure  4.3a. 
Upon  study,  this  figure  shows  the  densities  generated  by  the  VFVC  more  closely  fol¬ 
low  the  CTH  results  produced  with  the  Sesame  EOS.  This  is  contrary  to  the  results 
seen  in  Section  4.1.2,  Figure  4.1a  in  which  the  VFVC  generated  results  follow  the 
results  produced  in  CTH  using  the  Mie-Gruneisen  EOS.  This  discrepancy  is  a  result 
of  the  Mie-Gruneisen  EOS  coefficients,  which  were  the  same  for  the  VFVC  and  CTH 
using  the  Mie-Gruneisen  EOS  for  the  simulation  of  aluminum  impact  scenarios.  For 
iron,  however,  CTH’s  internal  e  phase  Mie-Gruneisen  EOS  coefficients  are  used  in  the 


Density  Behind  the  Compression  Wave 


Internal  Energy  Behind  the  Compression  Wave 


a.  Density  b.  Internal  Energy 


Pressure  Behind  the  Compression  Wave 


c.  Pressure 


Figure  4.3:  Comparison  of  Density,  Internal  Energy  and  Pressure  generated  by  the 
Finite  Volume  Code  and  CTH 

PTRAN  EOS,  while  the  VFVC  requires  the  input  of  the  e  Mie-Gruneisen  coefficients. 
In  the  case  of  iron,  the  coefficients  input  into  the  VFVC  happen  to  match  match  the 
material  properties  in  the  Sesame  EOS  better  than  they  do  for  the  PRTAN  EOS. 
This  discrepancy  highlights  the  need  for  to  validate  the  VFVC  with  analytical  data, 
experimental  data,  or  numerical  data  produced  by  a  reliable  hydrocode.  Figure  4.3b 
shows  the  numerical  values  of  internal  energy  generated  by  the  different  codes  and 
EOS,  once  again  are  all  quite  similar.  However  in  4.3b  it  is  apparent  that  at  impact 
velocities  grater  than  2.5  the  VFVC  generates  pressures  significantly  lower  than 
those  generated  by  CTH  using  both  EOS. 
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Table  4.5  contains  the  SD  comparing  the  numerical  values  produced  by  the 
VFVC,  and  CTH  using  the  Sesame  and  PTRAN  EOS.  The  SDs  of  0.0210  and  0.0268 
between  the  VFVC  generated  densities  and  those  generated  in  CTH  using  the  Sesame 
EOS  and  PTRAN  EOS  respectively,  are  less  than  one  percent  of  the  values  of  density 
being  generated  making  this  little  difference  acceptable.  The  SDs  of  the  differences 
in  internal  energies  generated  between  the  VFVC  and  CTH  using  the  Sesame  and 
PTRAN  EOS  are  18.285  and  17.8018  respectively.  These  values  of  SDs  maintain 
percentages  on  the  order  of  101  for  most  of  the  internal  energy  values  generated.  As 
seen  in  Figure  4.3b  it  is  at  the  higher  impact  velocities  where  the  difference  between  the 
VFVC  and  the  CTH  generated  values  of  internal  energy  have  the  greatest  differences 
in  value.  So  depending  on  the  level  of  accuracy  necessary,  there  exists  an  impact 
velocity  that  will  serve  as  an  upper  limit  for  the  applicability  of  the  VFVC  when 
modeling  impact  with  iron,  in  the  VFVC’s  present  form.  The  need  for  this  upper 
limit  becomes  quite  evident  when  considering  the  SD  for  the  values  of  the  numerically 
generated  pressures  are  1.5303  and  1.7815  GPa  for  the  difference  between  the  VFVC 
and  CTH  values  using  the  Sesame  and  PTRAN  EOS  respectively.  These  SD  values 
are  only  an  order  of  magnitude  smaller  than  the  values  of  the  pressures  produced  for 
most  of  the  range  of  impact  velocities,  which  is  a  good  percentage  of  the  generated 
values.  Looking  at  Figure  4.3b  it  is  apparent  that,  similar  to  the  internal  energy 
results,  the  differences  between  the  VFVC  and  CTH  generated  values  occur  at  the 
higher  impact  velocities. 


Table  4.5:  Comparison  of  the  Standard  Deviations  Between  the  VFVC,  Sesame 

EOS  and  Mie-Gruneisen  EOS  Generated  Values _ 


State 

Variable 

SD  of  VFVC  &  SD  of  VFVC  &  SD  of  Sesame  & 
Sesame  EOS  PTRAN  EOS  PTRAN  EOS 

Pressure  (GPa) 
Density  (%) 

Internal 

Energy 

1.5303  1.7815  0.3332 

0.0210  0.0268  0.0201 

18.267  17.802  4.2016 
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To  determine  if  the  difference  in  pressure  at  the  high  velocities  is  a  result  of 
the  Mie-Gruneisen  coefficients  used,  or  the  codes  inability  to  handle  multiple  phases, 
the  CTH  values  for  the  a  and  e  Mie-Gruneisen  coefficients  for  iron  were  used  in  the 
VFVC  for  impact  velocities  of  3500  ™  and  4000  ™ .  The  answers  generated  are  given 
in  Table  4.6.  From  Table  4.6  it  is  apparent  that  the  values  generated  by  the  VFVC 
with  the  Mie-Gruneisen  EOS  coefficients  from  CTH  are  2.6%  and  3.6%  lower  than  the 
corresponding  values  generated  by  CTH  for  the  impact  velocities  of  3500  and  4000 
This  indicates  the  differences  in  the  pressures  for  generated  by  VFVC  for  the  impact 
of  iron  is  caused  by  its  inability  to  handle  multiple  material  phases,  and  not  by  the 
Mie-Gruneisen  coefficients  used. 


Table  4.6:  Comparison  of  Pressures  Generated  at  Velocities  of  3500  —  and  4000  — 

1  S  S 

by  CTH  using  the  Sesame  and  Mie-Gruneisen  EOS  and  the  VFVC  Using  the  Mie- 
Gruneisen  Coefficients  from  Table  4.3  and  CTH _ 


Impact 

Velocity 

Sesame  EOS 

PTRAN  EOS 

VFVC 

Origional 

VFVC 

w/CTH  values 

3500  ^ 

s 

92  GPa 

92.4  GPa 

87.9  GPa 

89.6  GPa 

4000  m 

s 

112.1  GPa 

113.6  GPa 

106  GPa 

108  GPa 

Figure  4.4  shows  the  wave  profile  resulting  from  an  iron  on  iron  impact  with  an 
impact  velocity  of  612.7  "  at  a  time  of  1  /is.  Again,  the  added  dissipation  from  the 
Lax-Friedrich  flux  routine  is  smearing  out  the  elastic  precursor  in  the  results  generated 
by  the  VFVC  in  Figure  4.4c.  In  Figure  4.4b,  the  compression  waves  have  reached  the 
outside  boundaries  an  have  began  to  reflect  as  rarefaction  waves.  This  is  evident 
while  looking  at  two  things;  the  first,  is  the  sharp  cutoff  of  the  elastic  precursor  at 
the  boundaries,  second,  is  the  reduction  in  pressure  at  the  upper  ’’corners”  of  he 
compression  wave.  This  lower  pressure  is  due  to  the  beginning  of  the  rarefaction  wave 
relieving  the  pressure  caused  by  the  impact. 

Overall,  the  VFVC  appears  to  give  answers  quite  comparable  to  those  gener¬ 
ated  by  CTH,  especially  for  a  non-polymorphic  material.  The  necessity  of  validating 
the  initial  results  generated  by  the  VFVC  to  determine  if  the  proper  values  of  Mie- 
Gruneisen  coefficients  are  being  used,  has  also  been  established.  The  VFVC  code 
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Figure  4.4:  Comparison  of  Pressure  waves  generated  by  the  Finite  Volume  Code 

and  CTH 
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will  simulate  impact  with  polymorphic  material  better  if  a  subroutine  capable  of  han¬ 
dling  the  phase  transition  is  added.  Another  desired  addition  to  the  VFVC  is  an 
improvement  to  the  flux  routine. 
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V.  Comparison  of  CTH  and  NET  Hydrocode 

This  part  of  the  investigation  consists  of  comparing  the  numerical  results  generated 
in  by  Lu  and  Hanagud  in  Reference  [26]  to  the  results  generated  by  CTH  using  the 
Sesame  and  both  PTRAN  EOS.  The  impact  scenario  being  modeled  consists  of  an 
iron  target  1  mm  long  with  a  stress  load  of  50  GPa  applied  on  one  side  to  initiate  a 
stress  wave. 

In  their  simulations,  Lu  and  Hanagud  are  incorporate  a  slip  system  into  their 
grid  using  slip  angles,  </>,  of  0°  and  45°  degrees.  The  slip  system  of  0°  is  parallel  to  the 
axis  CTH  uses  when  modeling  uniaxial  impacts,  and  thus  is  automatically  accounted 
for  by  CTH.  To  introduce  a  slip  system  at  an  angle  of  45°  the  plane  strain  mesh  in 
Figure  5.1  is  utilized.  This  mesh  consists  of  a  layer  of  PMMA  inserted  into  the  target 
material  at  a  45°  angle.  The  rigid  boundary  layers  located  at  the  top  and  bottom  of 
the  mesh  prohibit  mass  from  flowing  out  of  the  mesh  limiting  any  strain  from  occurring 
in  the  vertical  direction,  effectively  creating  a  uniaxial  strain  scenario.  The  layer  of 
PMMA  acts  as  a  dampener  producing  lower  stresses  in  the  back  of  the  target  than  at 
the  front.  This  requires  a  higher  impact  velocity  to  account  for  the  dampening  effect 
of  the  PMMA.  In  the  plane  strain  mesh,  each  cell  measures  0.002  mm  by  0.002  mm, 
this  is  smaller  than  the  0.002  cm  cell  size  suggested  by  Szmerekovsky  [35],  however  it 
is  used  to  match  the  cell  size  used  by  Lu  and  Hanagud.  The  0.002  mm  cell  width  is 
also  used  in  the  uniaxial  case.  Overall,  the  uniaxial  mesh  is  2  mm  long  with  a  1  mm 
long  target  and  a  1  mm  long  flyer  plate.  The  length  of  the  grid  is  1000  cells  long. 
The  plane  strain  mesh  shares  the  same  overall  horizontal  dimensions  as  the  uniaxial 
grid  and  is  0.2  mm  tall,  the  overall  mesh  dimensions  are  1000  cells  in  the  horizontal 
direction  by  100  cells  in  the  vertical  direction. 

Lu  and  Hanagud  apply  stress  to  initialize  their  impact  scenario,  CTH  however, 
requires  a  specified  impact  velocity.  This  requires  the  determination  of  the  impact 
velocity  needed  to  generate  50  GPa  of  pressure.  To  do  this,  the  curves  in  Figure  3.14d 
are  fit  with  a  6th  degree  polynomial,  which  is  then  used  to  solve  for  an  initial  impact 
velocity  equal  to  2255.0  A  6th  degree  polynomial  is  used  because  it  came  closest 
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to  fitting  the  curve  in  the  region  of  interest.  This  velocity  initially  produces  a  stress 
smaller  than  the  requisite  50  GPa,  in  both  the  uniaxial  and  plane  strain  cases.  To 
fix  this  the  velocity  is  then  increased  incrementally  until  the  correct  impact  velocities 
are  determined.  The  final  impact  velocities  are  2265.05  ™  for  the  uniaxial  case  and 
2280  ™  for  the  plane  strain  case. 

The  1  mm  thick  flyer  plate  is  used  to  initiate  the  stress  wave  in  the  target.  To 
avoid  interference  from  wave  reflections  at  the  boundaries,  a  semi- infinite  boundary 
condition  is  used  at  the  far  left  and  right  boundaries. 

CTH  is  used  to  reproduce  Lu  and  Hanagud’s  results,  using  both  the  reversible 
and  irreversible  PTRAN  and  Sesame  EOS,  for  the  uniaxial  and  plane  strain  cases. 
The  results  generated  by  CTH  are  then  compared  to  those  generated  by  Lu  and 
Hanagud. 

5.1  Comparison  of  CTH  to  A  NET  Hydrocode 

In  this  section,  the  results  generated  by  Lu  and  Hanagud’s  NET  code  are  com¬ 
pared  to  the  results  generated  by  CTH  using  the  Sesame  and  both  PTRAN  EOS. 
Then,  any  similarities  or  differences  between  the  results  generated  by  the  two  hy¬ 
drocodes  results  will  be  analyzed. 

5.1.1  Uniaxial  Impact  with  oq  =  50 GPa  and  0  =  0°.  The  first  set  of  results 
to  be  compared  are  the  uniaxial  impact  results  where  oq  =  50 GPa  and  0  =  0°.  Figure 

5.2  contains  the  CTH  results  using  both  PTRAN  EOS  and  the  Sesame  EOS  and  the 
results  generated  by  Lu  and  Hanagud.  The  times  the  wave  profiles  are  captured 

Flyer  Plate  -*j  i  mm  Target 


Impact  Velocity  - ► 

1mm 

2mm 

\  1mm 

45°  Angle 

Figure  5.1:  Schematic  of  Mesh  with  45  Degree  Slip  System 
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at  are  1.1318  *  10“7,  1.1301  *  10“7,  1.13097  *  10“7  and  1.131  *  10“7  s  for  Lu  and 
Hanagud’s  results  and  those  generated  by  CTH  using  the  Sesame  and  the  reversible 
and  irreversible  PTRAN  EOS,  respectively.  Lu  and  Hanagud  use  y  as  the  horizontal 
axis  variable,  explaining  the  difference  in  notation  in  these  results  ang  those  presented 
for  the  plane  strain  case. 

Examining  the  results  contained  in  Figure  5.2,  small  differences  in  the  results 
generated  using  Lu  and  Hanagud’s  code  and  those  generated  by  CTH  using  the 
PTRAN  and  Sesame  EOS  may  be  ascertained.  First  the  difference  in  the  results 
generated  using  the  different  codes  are  presented.  The  most  obvious  difference  in 
results  generated  by  the  two  codes  is  the  amount  of  dispersion  error  (oscillations) 
in  Lu  and  Hanagud’s  results.  This  dispersion  error  is  caused  by  the  magnitude  of 
the  shockwave  developed  by  their  code  and  may  be  fixed  by  increasing  the  artificial 
viscosity  used  in  their  code.  Another  difference  in  the  results  generated  by  the  two 
codes  is  the  distance  between  the  beginning  of  the  elastic  precursor  and  the  begin¬ 
ning  of  the  plastic  wave,  of  which  CTH  has  the  larger  distance.  This  is  caused  by 
the  difference  in  plastic  wave  speed  generated  by  the  different  codes,  made  apparent 
by  the  extra  distance  that  the  plastic  wave  has  covered.  The  cause  of  the  different 
plastic  wave  speed  is  most  likely  the  use  of  different  material  parameters,  such  as  the 
reference  density  (po)  of  e  iron,  used  by  Lu  and  Hanagud  compared  to  those  contained 
in  CTH.  The  different  e  iron  material  properties  have  more  of  an  effect  than  the  a 
iron  properties  because,  as  is  evident  in  Figure  5.2,  the  e  iron  stress  wave  overtakes 
a  iron  wave  due  to  the  high  value  of  the  stress  applied.  The  only  similarity  in  the 
results  generated  by  the  different  codes  is  the  value  of  the  equilibrium  stress  behind 
the  compression  wave. 

The  results  generated  by  CTH  using  the  different  EOS  are  fairly  close.  The 
stress  generated  using  both  PTRAN  EOS  is  slightly  larger  than  50  GPa,  but  not 
significantly. 
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Unfortunately,  not  much  can  be  ascertained  from  the  results  generated  by  the 
different  hydrocodes  that  may  be  used  in  determination  of  the  different  codes  respec¬ 
tive  abilities  to  handle  NET. 

5.1.2  Uniaxial  Impact  with  <7o  =  50 GPa  and  (p  =  45°.  The  second  set 
of  results  to  compared  are  the  for  the  plane  strain  case  where  cy,  =  50 GPa  and 
<p  =  45°.  The  wave  reflections  caused  by  inserting  a  layer  of  PMMA  contribute  to  the 
generation  of  very  busy  wave  profiles.  This  necessitate  the  CTH  results  be  presenting 
using  wave  histories  instead  of  profiles.  A  histogram  created  by  placing  a  tracer 
point  just  behind  the  layer  of  PMMA  presents  the  wave  information  in  a  manner  in 
which  the  different  waves  generated  may  more  easily  be  distinguished  from  each  other. 
Figure  5.3  contains  the  CTH  results  using  the  Sesame  and  both  PTRAN  EOS  and 
the  results  generated  by  Lu  and  Hanagud.  Unlike  the  previous  impact  scenario,  the 
stress  waves  generated  in  this  case  consist  of  two  distinguishable  waves.  As  Lu  and 
Hanagud  explain  the  initial  wave  is  a  elastic  wave  produced  by  the  e  iron,  a  plastic 
wave  then  follows  behind  the  elastic  wave  [26] .  The  different  waves  are  labelled  in  the 
graph  generated  by  Lu  and  Hanagud  in  Figure  5.3c. 

The  histograms  produced  by  CTH  show  several  waves.  The  first  of  these  waves 
have  a  peak  stress  between  10  GPa  and  15  GPa.  This  wave  begins  at  the  PMMA/iron 
boundary  and  propagates  outward  beginning  at  he  time  of  impact.  It  is  a  result  of  the 
method  used  to  create  the  45°  slip  system  in  the  CTH  model,  and  is  not  a  real  solution, 
thus  it  can  be  ignored.  The  second  wave  corresponds  to  the  elastic  wave  in  Lu  and 
Hanagud’s  results.  The  third  wave  is  the  plastic  wave  covered  with  oscillations. 

The  similarities  in  the  compression  waves  generated  using  Lu  and  Hanagud’s 
code  and  CTH  include  the  production  of  separate  elastic  and  plastic  waves  and  the 
production  a  a  50  GPa  equilibrium  stress  value  behind  the  plastic  wave.  CTH’s 
production  of  elastic  and  plastic  waves  validates  the  simulation  of  a  slip  system  by 
inserting  a  layer  of  softer  material  in  the  PMMA  at  the  desired  slip  angle.  The 
production  of  multivalent  equilibrium  stresses  behind  the  plastic  waves  is  expected, 
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because  the  equilibrium  thermodynamics  included  in  both  codes  should  cause  them 
to  generate  the  same  equilibrium  value  of  stress. 

Another  similarity  observable  in  Figures  5.3a  and  5.3b  are  the  stress  wave  his¬ 
tories  produced  using  the  Sesame  and  both  PTRAN  EOS.  This  indicates  any  ability 
CTH  may  have  in  modeling  the  NET  of  this  scenario  are  not  generated  by  the  EOS 
and  are  most  likely  the  product  of  CTH’s  ability  apply  constitutive  relationships. 

Several  differences  in  the  answers  generated  using  the  different  codes  are  ob¬ 
served  in  Figure  5.3.  These  differences  include  the  shape  and  magnitude  of  the  elastic 
wave  and  oscillations  evident  in  CTH’s  representation  of  the  plastic  wave.  The  large 
oscillations  are  likely  caused  by  excessive  interference  from  wave  reflections  resulting 
from  inserting  a  layer  of  PMMA  in  the  iron  to  produce  the  slip  system. 

Comparing  the  difference  in  magnitude  and  shape  of  the  elastic  wave  allows 
for  the  comparison  of  the  different  hydrocodes’  abilities  to  model  the  NET  of  an 
impact  system.  In  Figure  5.3  the  elastic  wave  in  Lu  and  Hanagun’s  results  reaches  a 
peak  stress  of  approximately  46  GPa  before  sharply  relaxing  to  a  stress  of  just  under 
30  GPa.  This  contrasts  with  the  elastic  wave  generated  by  CTH  peaking  at  a  stress 
of  approximately  37  GPa  before  relaxing  to  a  stress  of  30  GPa.  Lu  and  Hanagud’s 
hydrocode  also  has  a  faster  relaxation  time  than  seen  in  the  CTH  results,  indicative 
of  their  code’s  ability  to  handle  the  non-equilibrium  region  immediately  behind  a 
shockwave  where  NET  conditions  abound. 


The  difference  in  peak  stress  value  is  a  result  of  Lu  and  Hanaguds  hydrocode 
to  account  for  a  non-equilibrium  stress  component  as  seen  in  Equation  (2.93)  shown 
here. 


In  Equation  (2.93)  the  total  stress  tensor  is  composed  of  an  equilibrium  stress  and  a 
non-equilibrium  stress,  a  total  =  &eq  +  &neq .  In  this  case,  the  9  GPa  difference  in  the 
peak  stresses  is  produced  by  the  non-equilibrium  stress  component,  aneq  included  in 
Lu  and  Hanagud’s  hydrocode.  This  stress  component  is  a  found  using  the  differential 
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equation  given  in  Equation  (2.89),  shown  here  for  convenience. 


+  Vm  V% 

that  includes  the  strain  rate,  V£t,  and  the  relaxation  time  of  the  viscous  processes, 
ru  [26].  These  results  suggest  the  Lu  and  Hanafud’s  hydrocode  is  capable  of  modeling 
system  in  non-equilibrium  states. 

While  not  accounting  directly  for  iy,  CTH  does  use  strain  rate  in  the  Zerilli- 
Armstrong  constitutive  relationship,  allowing  it  to  generate  the  first  stress  wave  and  to 
account  for  some  relaxation  of  the  material  after  the  stress  wave  passes.  This  is  shown 
in  Figures  5.2a  and  5.2b.  This  indicates,  while  CTH  does  not  use  the  irreversible, 
non-equilibrium  variables  from  the  theories  of  IVT  and  EIT,  it  is  still  capable  of 
accounting  for  the  irreversible  effects  of  material  properties  including  stresses  and 
strains.  In  addition,  CTH  is  also  able  to  account  for  the  irreversibilities  of  phase 
changes. 
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VI.  Conclusions 


6.1  Conclusions 

This  work  had  three  separate  objectives.  The  first  objective  was  to  deter¬ 
mine  the  similarities  and  differences  between  the  tabular  Sesame  EOS,  and  the  semi- 
empirical  Mie-Gruneisen  and  PTRAN  EOS.  This  objective  provides  guidance  on  the 
application  of  the  EOS  to  impact  scenarios  using  both  CTH  and  the  FVC.  The  second 
objective  included  validating  the  FVC  by  comparing  its  answers  to  those  generated  in 
CTH  using  the  Sesame  and  Mie-Gruneisen  EOS.  This  validation  is  a  key  step  in  the 
development  of  the  FVC,  which  has  the  possibility  being  developed  into  a  hydrocode 
capable  of  modeling  non-equilibrium  thermodynamics.  The  final  objective  involved 
comparing  CTH  results  to  those  generated  by  a  program  capable  of  handling  NET  to 
determine  CTH’s  ability  to  model  NET. 

Several  conclusions  were  made  in  the  comparison  of  the  Sesame  and  Mie-Gruneisen 
EOS.  The  first  and  most  important  conclusion  for  the  further  production  of  the  VFVC, 
stem  from  the  similar  values  of  internal  energy  and  pressure  generated  by  CTH  us¬ 
ing  either  the  Sesame  or  Mie-Gruneisen  and  PTRAN  EOS.  These  results  lead  to  the 
conclusion  the  Mie-Gruneisen  EOS  produces  the  similar  results  to  those  generated  by 
the  Sesame  EOS  for  the  same  impact  velocities.  However,  the  Mie-Gruneisen  EOS 
cannot  produce  the  same  results  as  the  Sesame  EOS  for  any  value  of  density,  or  at 
temperatures  above  the  Te  in  the  Sesame  EOS  for  the  material  being  modeled.  It  was 
also  shown  that  all  the  EOS  are  more  sensitive  to  density  then  to  temperature.  An¬ 
other  finding  is,  the  Sesame  and  Mie-Gruneisen  EOS  will  produce  different  densities, 
while  producing  the  similar  pressures.  This  difference  may  be  corrected  by  changing 
the  Mie-Gruneisen  coefficients  used,  however  this  will  affect  the  pressures  generated. 
This  highlights  the  need  to  validate  numerically  generated  answers  using  experimen¬ 
tal  data  or  analytical  solutions,  to  assure  generated  values  for  a  parameter  of  interest 
are  correct.  Finally,  and  most  importantly  in  the  determination  of  CTH’s  ability  to 
model  NET,  the  Sesame  EOS  is  shown  to  generate  irreversibilities  in  regions  of  phase 
changes. 
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The  major  conclusion  stemming  from  the  second  objective  is  the  validation 
of  the  FVC  using  the  hydrocode  CTH.  The  FVC  produced  values  of  pressure  and 
internal  energy  similar  to  the  numerical  answers  generated  by  CTH  using  the  Sesame 
and  Mie-Gruneisen  EOS  at  the  same  impact  velocities  in  a  non-polymorphic  material. 
For  a  polymorphic  material,  however,  the  pressures  and  internal  energies  generated  by 
the  FVC  were  considerably  different  at  high  speeds.  This  indicates  the  need  to  add  a 
multi-phase  capability  to  the  FVC  in  order  for  the  code  to  accurately  model  impacts 
involving  polymorphic  materials.  The  need  for  an  flux  routine  that  can  accurately 
produce  the  HEL  caused  by  impacts  was  also  demonstrated. 

In  the  comparison  of  results  generated  by  CTH  to  those  generated  by  Lu 
and  Hanagud’s  NET  hydrocode  two  conclusions  were  reached.  The  first,  Lu  and 
Hanagud’s  code  is  able  to  model  a  system  outside  an  equilibrium  state  using  irre¬ 
versible,  non-equilibrium  thermodynamics  from  the  taken  from  theories  of  IVT  and 
EIT.  The  second,  and  maybe  more  significant,  is  CTH  is  able  to  account  for  some 
of  the  same  irreversible  parameters  including  stress,  strain  and  strain  rate  contained 
in  the  constitutive  relationship  and  the  phase  changes  in  polymorphic  materials  con¬ 
tained  in  the  EOS. 

6.2  Recommendations 

In  this  investigation  the  ability  of  the  Mie-Gruneisen  and  Sesame  EOS  to  model 
uniaxial  impact  scenarios  was  investigated.  There  are  more  EOS  to  be  investigated 
for  the  possible  use  in  a  hydrocode  that  can  handle  NET.  Another  investigation  not 
addressed  in  this  work,  is  to  comparison  of  EOS  for  alloys  instead  of  those  for  pure 
metals. 

The  study  of  NET  is  currently  of  great  interest  and  is  heavily  reported  on.  To 
further  determine  CTH’s  ability  to  model  NET,  more  comparisons  should  be  made 
with  the  numerical  results  presented  in  journals  and  conferences. 
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Recommendations  for  the  further  development  of  the  FVC  include  alterations 
to  the  flux  schemes,  validation  of  the  FVC’s  constitutive  model  and  adding  an  ability 
to  handle  multiple  phases.  For  the  flux  scheme  change  it  is  recommended  the  Lax- 
Friedrich’s  flux  scheme  be  abandoned,  then  using  a  different  1st  order  flux  scheme. 
Another  option  is  including  higher  order  (3rd  and  4th  order)  flux  schemes,  for  example 
the  4th  order  essentially  non-oscillatory  (ENO)  scheme  detailed  in  Udaykumar  et 
al.  [38].  The  constitutive  model  should  be  validated  allow  for  accurate  plastic  flow 
results.  This  is  a  key  step  in  the  development  of  a  code  that  can  model  NET,  because 
the  plastic  deformation  is  a  source  of  irreversible  dissipation.  Also,  if  impacts  of 
polymorphic  materials  are  going  to  be  modeled,  the  code  needs  to  be  adopted  to 
handle  multi-phase  materials. 
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Appendix  A.  Conversion  of  CFD  Code  into  a  Hydrocode 

This  Appendix  contains  portions  of  a  document  written  describing  the  equations 
and  procedure  applied  to  an  one  dimensional  CFD  code  intended  to  model  the  fluid 
dynamics  in  a  shock  tube  for  the  purpose  of  transforming  it  into  an  one  dimensional 
hydrocode  capable  of  modeling  uniaxial  impacts.  The  document  is  intended  to  be  a 
stand  alone  document  with  its  own  bibliography,  albeit  a  small  one.  The  document 
also  assumes  the  reader  is  familiar  with  computational  mechanics.  However,  after 
reading  this  thesis,  it  should  not  be  difficult  to  understand. 
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Conversion  of  1-D  Shock  Tube  Code  to  a  1-D  Hydrocode 


I.  Introduction 


This  paper  describes  the  conversion  of  a  1-D  Shock  Tube  CFD  code  into  a  1-D 
Hydrocode  capable  of  modeling  the  propagating  waves  caused  by  the  impact  of  two 
metal  bars  of  the  same  material.  In  addition  to  the  conservation  equations  already 
contained  in  the  shock  tube  code,  constitutive  equations  accounting  for  the  strength  of  the 
material  will  be  included  as  well  as  a  source  term.  In  effect,  the  hyperbolic  partial 
differential  equation  being  solved  is  : 


dQ  dF_ 
dt  dx 


=  5(0 


In  addition  to  the  constitutive  equations  and  source  terms,  an  equation  of  state  will 
be  utilized  to  relate  the  internal  energy  and  density  of  the  material  to  its  pressure. 


II.  Theory  and  Equations 


GOVERNING  EQUATIONS 

The  governing  equations  used  in  the  new  code  consist  of  the  conservation  laws  of 
mass,  momentum  and  energy  (equations  1,  2  and  3),  and  constitutive  equations  that 
account  for  the  strength  of  the  material.  Equation  4  accounts  for  the  plastic  stain  seen  by 
the  material,  while  equation  5  accounts  for  the  deviatoric,  or  shear  stress  seen  by  the 
material  . 


^+5(a)=0 

dt  dx 
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Though  equation  4  is  completely  uncoupled  from  the  rest  of  the  governing  equations,  the 
strain  experienced  in  an  impact  is  of  interest,  thus  explaining  its  inclusion  in  this 
assignment.  The  vector  of  conserved  variables  and  the  flux  vector  are  determined  from 
the  governing  equations  and  are2: 
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EQUATIONS  OF  STATE 


As  previously  mentioned  the  Mie-Gruneisen  EOS  is  used  to  relate  pressure  to  the 
internal  energy  (e)  and  the  specific  volume  (V  =  1/  p).  The  form  of  the  Mie-Gruneisen 
equation  used  in  this  code  is2: 


p(e,V) 


c02(y0-u)  j<y) 
[V0-s(V0-V)f  v 


2 


c0(Vo-V) 

Uo-j(V0-V)J 


Where  in  copper  Co  is  the  bulk  speed  of  sound  at  zero  pressure,  s  is  the  slope  of  the 
Hugoniot1  and  T  is  the  Gruneisen  parameter.  As  the  waves  flow  through  the  materials  T 
is  dependent  on  the  density  of  the  material  and  is  subject  to  the  following  relationship  : 

r=v( &}  -  r°Po 

V  Jy  P 

These  parameters  are  also  used  to  calculate  the  wave  speed  of  the  plastic  wave  moving 
through  the  material  following  the  equation  : 

Up  =c0  +  sup 

where  up  denotes  the  particle  velocity.  To  determine  speed  of  the  elastic  wave,  the 
*  ^ 
following  relationship  is  used  : 


where  K  is  the  bulk  modulus  and  is  equal  to5: 

K  =  pc02 

When  using  the  Mie-Gruneisen  EOS,  the  speed  of  sound  in  the  metal  is  found  using  the 
relationship2: 
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c 
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P  +  Poco 


P+(s-r0)(p-p0) 

[p-s(p-P0)J 


The  value  of  the  total  Cauchy  stress  tensor  produced  by  the  elastic  and  plastic 
waves  is  determined  with  the  realationship2: 


<t  —  r  +  pS 

lj  lj  ^  lj 


TIME  INTEGRATION 

In  order  to  integrate  the  PDE  in  time,  the  it  is  rearranged  to  resemble  : 


where 


L(Q)  =  -^-^  +  D(Q) 

Ax 

and  D  is  the  operator  designating  the  source  terms.  The  source  term  is  added  to  the 
residuals  before  integrating  through  time. 

Due  to  the  small  time  steps  made  necessary  by  the  high  speed  of  sound  in  a  metal, 
three  different  methods  of  time  integration  were  utilized  for  this  report,  the  Euler  Explicit 
method  used  in  the  previously  unaltered  code  and  two  Runge  Kutta  Schemes.  The  first 
Runge-Kutta  scheme  is  a  three  stage,  3rd  order  TVD  Runge-Kuttal  Scheme,  optimized  by 
Shu  and  Osher  for  handling  ENO  shock  capturing  schemes  and  applied  by  Udaykumar 
et  al.  to  wave  propagation  through  materials.  Though  an  ENO  scheme  is  not  being 
utilized  in  this  code,  the  decision  was  made  to  include  it  to  determine  if  it  has  any 
advantage  in  solving  this  type  of  problem,  thereby  improving  the  numerical  solutions 
generated.  This  Runge-Kutta  scheme  takes  the  following  form: 

M(1)=M(0)+ArL(n(0)) 

M(2)=2M(o)+IM(D+IAfL(M(i)) 

4  4  4  v  ' 

„®=V>+Am+-A/L(«ra) 

3  3  3  v  ' 

th  4 

A  four  stage,  4  order  Runge-Kutta  is  also  utilized  : 
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Md)  =m(0)+^(l(o)) 

m(2)=m(0)+^(l(1)) 

M(3)=u(0)+At(L(2)) 

w(n+1)  =  um  +  y(L(0)  +  2L(1)  +  2  L(2)  +  L(3) ) 

For  1st  order  spatial  integration  the  Lax-Friedrichs  routine  in  the  previously 
existing  code  was  utilized.  In  order  to  get  second  order  spatial  results  the  MUSCL 
scheme  with  Min  Mod  limiters  is  utilized4. 

III.  Implementation 

The  following  is  a  description  of  the  code  changes  implemented  how  some 
data  it  stored.  A  detailed  list  of  the  code  additions  and  alterations  is  located  in  Appendix 
A. 

The  implementation  of  the  code  consists  of  reading  in  and  initializing  the  material 
properties  of  copper,  updating  the  Q  vector  to  handle  the  additional  governing  equations 
and  material  properties,  updating  the  flux  routine,  incorporating  the  source  terms  into  the 
right  hand  side  module,  incorporating  the  new  Runge-Kutta  time  integration  techniques 
into  the  solve  module,  and  updating  the  output  module. 

The  material  properties  are  incorporated  into  to  the  input  deck,  along  with  a 
hydrodynamic  operator,  spatial  integration  designator  and  time  integration  technique 
designator.  These  designate  the  program  to  be  either  a  pure  hydrodynamic  code  or  a  code 
with  the  ability  to  handle  elastic-plastic  conditions  in  addition  to  specifying  the  spatial 
and  time  order  of  accuracy.  A  maximum  time  value  is  used  to  tell  the  program  when  to 
stop.  The  material  properties,  hydrodynamic  operator  and  max  time  are  read  in  inside  the 
initial  module  and  are  stored  in  the  runtime  module.  The  initial  module  is  also  used  to 
initialize  the  material  properties  and  the  initial  velocities  in  the  grid.  In  both  the 
hydrodynamic  and  elastic-plastic  settings  the  plastic  strain  and  deviatoric  stresses  are 
initialized  to  zero. 

In  addition  to  the  existing  governing  equations  already  contained  in  the  Q  vector 
(mass,  momentum  and  energy),  the  plastic  strain  (PS)  and  the  deviatoric  stress  (Sx)  have 
been  added.  The  values  of  the  wave  speeds  and  the  total  values  of  the  Cauchy  stress 
tensor  are  also  carried  in  the  Q  vector  for  convenience,  thought  if  computing  time  is  an 
issue,  these  values  may  also  be  calculated  at  the  end  of  the  code.  The  definitions  of  the 
new  Q  vectors  are  contained  in  the  constant  module,  and  changes  to  the  Q  vector  itself 
are  located  throughout  the  code. 

The  Lax-Friedrich  flux  routine  is  utilized  to  perform  the  spatial  integration. 
Modifications  to  the  Lax-Friedrichs  subroutine  include;  adding  the  additional  flux  terms, 
implementing  the  Mie-Gruneisen  Equation  of  state  to  solve  for  the  pressure,  and 
incorporating  the  speed  of  sound  through  the  metal.  Since  the  pressures  and  speed  of 
sound  for  each  cell  is  different,  these  values  were  treated  as  primitive  values  and 
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calculated  each  time  for  each  cell  every  time  step.  For  2nd  order  spatial  results  the 
MUSCL  subroutine  with  the  Min-Mod  limiter  is  utilized. 

The  source  terms  are  included  inside  the  RHS  module,  where  an  if  statement  is 
utilized  to  either  set  the  source  terms  to  zero  for  the  hydrodynamic  solver,  or  to  use  the 
previously  derived  source  term  values  in  the  case  of  the  elastic-plastic  solver.  The  values 
for  the  individual  cell  velocities  and  deviatoric  stress  tensors,  including  the  ghost  cells, 
are  calculated  in  an  initial  do  loop  and  are  stored  in  arrays.  The  gradients  and  source 
terms  are  then  calculated  in  a  separate  do  loop.  At  the  end  of  the  RHS  module,  the  source 
term  is  added  to  the  residual  vector,  the  sum  of  which  is  then  passed  into  the  solve 
module  where  the  time  integration  is  performed. 

The  new  Runge-Kutta  integration  techniques  are  incorporated  into  the  solve 
module.  Both  integration  techniques  are  coded  in  a  series  of  successive  stages  since 
neither  really  allow  themselves  to  utilize  a  repeated  do  loop.  For  this  reason,  a  right  hand 
side  array  is  allocated  to  store  the  residuals  for  each  stage,  to  be  used  specifically  by  the 
four  stage  Runge-Kutta.  An  allocated  array  is  also  used  to  store  the  original  values  of  the 
Q  vectors  (designated  Q  old).  The  cells  pressures,  wave  speeds,  and  total  Cauchy  stress 
values  are  also  calculated  in  the  solve  module  using  the  equations  listed  in  the  theory 
section  of  this  report.  As  in  the  original  code,  the  time  step  is  calculated  in  each  cell  with 
the  minimum  cell  time  step  designated  the  global  time  step.  The  calculation  of  the  speed 
of  sound  used  in  the  time  step  was  updated  to  calculate  the  speed  of  sound  for  a  wave 
moving  through  a  metal.  Then  due  to  the  high  magnitude  of  the  speed  of  sound  through  a 
metal,  the  global  time  stem  is  divided  by  100. 

The  last  major  code  changes  include  modifying  the  output  module  to  output 
additional  desired  information  to  a  data  (.dat)  file. 
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Appendix  B.  CTH  Input  Deck 


This  Appendix  contains  the  CTH  Input  deck  used  to  model  the  impact  of  an  iron 
target  with  a  45°  slip  system  introduced.  This  is  the  input  deck  used  to  generated 
the  results  in  Section  5.1. 


Ill 


*eor*  genin 

Flyer  Plate  Test  With  45  Degree  Slip  System 
control 
nuiip 
ep 

vpsave 

endcontrol 

mesh 

block  1  georn=2dr  type=e  *  2dr  is  two  dimensional  rectangle  plain  strain 
x0=-0.1 

xl  n=1000  w=0.2  dxf=0.0002 
endx 


y0=-0.01 

yl  n=100  w=0.02  dyf=0.0002 
endy 

endb 

endmesh 

insertion  of  material 
block  1 

package  projectile  *define  projectile 
material  1 
numsuh  50 

xvel  2280e2  *  change  only  the  first  number,  leave  'e2'  this  converts  m/s  to  cm/s 

insert  box 
pi  -0.1, -0.01 
p2  0,  0.01 

endinsert 

endpackage 

package  target  *  insert  target 
material  2 
numsub  50 
insert  box 
pi  0,  -0.01 
p2  0.04,  0.01 
endinsert 
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delete  TRIANGLE 
PI  0.02,  -0.01 
P2  0.04,  -0.01. 

P3  0.04,  0.01 
ENDdelete 

endpackage 

package  PMMA  *  Insert  Target 
material  3 
numsub  50 

INSERT  PARALLELOGRAM 
PI  0.02,  -0.01 
P2  0.04,  0.01 
P3  0.03,  -0.01 
ENDINSERT 
endpackage 
package  back 
material  2 
numsub  50 

insert  TRIANGLE 

PI  0.03,  -0.01 
P2  0.05,-0.01. 

P3  0.05,  0.01 
ENDinsert 
endpackage 
package  back 
material  2 
numsub  50 
insert  box 

pi  0.05,  -0.01 

p2  0.1,0.01 

endinsert 

endpackage 

endblock 

endinsertion 
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edit 
block  1 
expanded 
endblock 
endedit 

tracer  *  add  tracer  points 

add  0.01  0 
add  0.035  0 
add  0.045  0 
add  0.05  0 
add  0.055  0 
endt 

eos 

*  Information  for  materials 
MAT1  ses  iron 
MAT2  ses  iron 
MAT3  ses  PMMA_RP 
endeos 

epdata 
mix  3 

matep  1  ZE  iron 
matep  2  ZE  iron 
matep  3  JO  LEX  AN 
vpsave 
lstrain 
endep 
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*eor*  cthin 

Flyer  Plate  Test  With  45  Degree  Slip  System 
control 
mmp 

tstop  =  2.0e-7  *  you  may  need  to  increase  the  stop  time  to  reach  end  of  event 

endc 


Convct 

convection=l 

interface=high_resolution 

endc 


edit 

shortt 


tim  0.0, 
ends 

dt=  1.0 

longt 
tim  0.0, 
endl 

dt=  1.0 

plott 
tim  0.0 

dt  =  1.0e-8 

*  this  sets  the 

endp 

histt 

tim  0.0, 

dt  =  0.1e-8 

*  this  sets  the 

htracer  all 
endh 
ende 


time  step  it  plots  profile  data 
time  step  it  plots  history  data 


boundary 
bhy 
bl  1 

bxb  =  1  ,  bxt  =  1 
byb  =  0,  byt  =  0 
endb 
endh 
endb 


*  semi-infinite  boundaries  at  x  max  and  min 
*rigid  boundaries  at  y  max  and  min 


cellthermo 

mmp3  *  This  was  recommended  by  Eglin  and  appears 

*  to  give  good  results  as  well, 
ntbad  5000000 
endc 
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Appendix  C.  CTH  Input  Deck 

This  Appendix  contains  the  spreadsheets  used  to  generate  the  graphs  for  impact 
scenarios  comparing  the  Mie-Gruneisen  and  Sesame  EOS  using  aluminum,  the  com¬ 
parison  of  the  PTRAN  and  Sesame  EOS  using  iron  and  the  validation  of  the  VFVC. 
The  values  contained  in  these  spreadsheets  were  generated  using  both  CTH  and  the 
VFVC. 
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a.  Mie-Gruneisen  EOS  Generated  Results 


Impact 

m/s 

Stress 

Gpa 

Temp 

K 

Density 

kg\mA3 

Int  Energy 
kJ/kg 

Velocity 

m/s 

Pressure 

Gpa 

Strain 

10A  -3 

500 

4.04167 

332.188 

2.8775 

300.938 

250 

3.78571 

21.0714 

1000 

8.41667 

372.75 

3.00658 

395 

500 

8.16667 

44.1667 

1500 

13.2857 

425.625 

3.12692 

550 

742.875 

13 

66.25 

2000 

18.5 

497 

3.245 

765 

1000 

18.4 

86.6667 

2500 

24.25 

592.647 

3.35938 

1041.67 

1250 

24.125 

104.167 

3000 

30.3125 

718.782 

3.47083 

1377.78 

1500 

30.3125 

127.5 

3500 

37.2727 

878.125 

3.575 

1778.57 

1750 

37.0455 

145 

4000 

43.8889 

1079.17 

3.67647 

2240 

2000 

43.9474 

158.929 

4500 

51.875 

1322.22 

3.78571 

2750 

2250 

51.875 

175 

5000 

60 

1607.14 

3.89091 

3343.75 

2500 

60 

178.571 

5500 

68.0952 

1941.67 

3.98095 

3964.29 

2750 

68.0952 

166.667 

6000 

77.7 

2315.78 

4.07292 

4666.67 

3000 

77.5 

151.866 

b.  Sesame  EOS  Generated  Results 


Impact 

m/s 

Stress 

Gpa 

Temp 

K 

Density 

kg\mA3 

Int  Energy 
kJ/kg 

Velocity 

m/s 

Pressure 

Gpa 

Strain 

10A  -3 

500 

3.92308 

332.233 

2.82067 

283 

250 

3.67857 

21.1234 

1000 

8.17391 

371.923 

2.94211 

376.25 

500 

7.92 

44.16667 

1500 

12.8663 

423.333 

3.05962 

531.944 

760 

12.6875 

66.25 

2000 

18.2105 

493.684 

3.17105 

745.238 

1000 

17.9 

85 

2500 

23.75 

587.5 

3.28125 

1025 

1250 

23.625 

104.1667 

3000 

29.8438 

709.091 

3.38887 

1361.2 

1500 

29.8438 

122.5 

3500 

36.4583 

864.706 

3.49091 

1757.07 

1750 

36.4583 

138.235 

4000 

43.4783 

1050 

3.59 

2260 

2000 

43.5938 

157.143 

4500 

51.1111 

1280 

3.68824 

2726.16 

2250 

50.8333 

170.522 

5000 

58.75 

1540 

3.7875 

3312.5 

2500 

58.75 

181.25 

5500 

67.1429 

1841.67 

3.875 

3964.29 

2750 

67.1429 

183.75 

6000 

75.8333 

2181.82 

3.97037 

4659.09 

3000 

76,25 

157.83 

c.  VFVC  Generated  Results 


Impact 

m/s 

Stress 

Gpa 

Temp 

K 

Density 

kg\mA3 

Int  Energy 
kJ/kg 

Velocity 

m/s 

Pressure 

Gpa 

Strain 

10A  -3 

500 

3.922 

n/a 

2.887 

299.6 

249.9 

3.898 

3.001 

1000 

8.292 

n/a 

3.01 

393.5 

500 

8.266 

5.778 

1500 

13.12 

n/a 

3.13 

549.4 

750 

13.1 

8.384 

2000 

18.42 

n/a 

3.246 

768.4 

1000 

18.39 

10.84 

2500 

24.18 

n/a 

3.359 

1049 

1250 

24.15 

13.17 

3000 

30.4 

n/a 

3.469 

1329 

1500 

30.37 

15.35 

3500 

37.08 

n/a 

3.575 

1798 

1750 

37.05 

17.43 

4000 

44.23 

n/a 

3.678 

2266 

2000 

44.2 

19.5 

4500 

51.83 

n/a 

3.779 

2797 

2250 

51.8 

21.35 

5000 

59.91 

n/a 

3.876 

3390 

2500 

59.78 

23.12 

5500 

68.44 

n/a 

3.971 

4045 

2750 

68.4 

24.88 

6000 

77.43 

n/a 

4.064 

4766 

3000 

77.38 

26.47 

Figure  C.l:  Spreadsheets  Generated  Using  Aluminum  Impacts 
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a.  PTRAN  EOS  Generated  Results 


Impact 

m/s 

Stress 

Gpa 

Temp 

K 

Density 

kg\mA3 

Int  Energy 
kJ/kg 

Velocity 

m/s 

Pressure 

Gpa 

Strain 

400 

7.8333 

324.5 

8.1775 

152.717 

200 

7.579 

0.021667 

500 

9.875 

332.031 

8.25882 

164.063 

250 

9  625 

0  0272 

612.7 

12.2 

341.19 

8.35 

180 

306.25 

12.087 

0.033 

671.1 

13.0909 

344.25 

8.41154 

185.882 

328.571 

12.8182 

0.0345 

991.6 

17.2727 

391.667 

8.9142 

291.818 

500 

17.1591 

0.075 

997 

17.3636 

391.667 

8.91176 

292.727 

500 

17.2727 

0.074546 

1150 

20.5263 

421.333 

9.01429 

326.316 

577.778 

20.5263 

0.082727 

1247 

22.75 

440 

9  08667 

348  75 

626  087 

22.75 

0.086316 

1292 

23.8158 

452.005 

9.12308 

361.25 

643.75 

23.8158 

0.088 

1292 

23.75 

449.565 

9.12 

360.714 

646.667 

23.75 

0.088261 

1396 

26.25 

473 

9  19286 

390.278 

700 

26.25 

0.091765 

1557 

30.3125 

511.25 

9.30769 

440.625 

776.923 

30.3125 

0.097619 

1567 

30.6667 

517  105 

9.31304 

445.313 

785.714 

30.5 

0.096739 

1871 

38.8636 

623.143 

9.54 

565 

940 

38.8636 

0.1075 

1887 

39.2 

628.846 

9.54545 

570 

945 

39.2 

0.109091 

1900 

39.7727 

634.615 

9.56 

577.5 

960 

39.5455 

0.102273 

1908 

40 

638.462 

9.56364 

582.5 

955 

40 

0.117273 

2500 

57.5 

952.381 

9.9625 

891.67 

1250 

57.5 

0  12778 

3000 

74.2857 

1304.35 

10.2917 

1234.78 

1500 

74.2857 

0.145588 

3500 

92.35294 

1750 

10.6094 

1619.05 

1750 

92.35294 

0.151389 

4000 

111.7241 

2269.231 

10.8971 

2069.44 

2000 

113.6 

0.134615 

b.  Sesame  EOS  Generated  Results 


Impact 

m/s 

Stress 

Gpa 

Temp 

K 

Density 

kg\mA3 

Int  Energy 
kJ/kg 

Velocity 

m/s 

Pressure 

Gpa 

Strain 

400 

7.88 

325.395 

8.2026 

152.935 

200 

7.579 

0.021786 

500 

9.875 

333  167 

8.28333 

164.412 

250 

9.6875 

0  0272 

612.7 

12.3636 

343 

8.3731 

180.278 

306.6738 

12.09091 

0.303 

671.1 

13.1818 

336.667 

8.41154 

187.059 

328.571 

12.9091 

0.035556 

991.6 

16.9 

363.28 

8.9871 

299.091 

495.24 

16.8 

0.090625 

997 

17.0476 

363.077 

8  9889 

300 

496.875 

17.0476 

0.091177 

1150 

20.05556 

392.5 

9  08333 

329.412 

575 

20.14706 

0.091667 

1247 

22.2619 

415.714 

9.14925 

352  632 

623.8095 

22.25 

0.091875 

1292 

23.2813 

425 

9.18333 

364.474 

640 

23.4375 

0.094118 

1292 

23.3021 

423.077 

9.17761 

363.158 

646.154 

23.4211 

0.091579 

1396 

25.8333 

447.619 

9.25385 

391.667 

700 

25.8333 

0.09375 

1557 

29.8214 

485.882 

9.37143 

441.071 

778.261 

29.8333 

0.095263 

1567 

30.19231 

490.588 

9.38 

444.643 

775 

30 

0.095238 

1871 

38.3333 

583.929 

9.6 

560.526 

936  842 

38.3333 

0.1 

1887 

38.6957 

590 

9.61111 

566.667 

938.889 

38.6957 

0.10087 

1900 

39.0476 

598.214 

9.62222 

573.684 

947.368 

39.0476 

0.101191 

1908 

39  3182 

600 

9.63 

576.316 

952.632 

39  1667 

0.101667 

2500 

57  222 

897.368 

10.0227 

897.12 

1250 

57  222 

0.125 

3000 

73.75 

1208 

10.3594 

1227.27 

1500 

73.75 

0.1417 

3500 

92.5 

1590.91 

10.6563 

1625 

1750 

92 

0.15333 

4000 

112.308 

2018.52 

10.9688 

2073.529 

2000 

112.143 

0.13875 

c.  VFVC  Generated  Results 


Impact 

m/s 

Stress 

Gpa 

Temp 

K 

Density 

kg\mA3 

Int  Energy 
kJ/kg 

Velocity 

m/s 

Pressure 

Gpa 

Strain 

400 

7.733 

n/a 

8.205 

153.4 

200 

7.683 

0.0136 

500 

9.809 

n/a 

8278 

163.6 

250 

9.75 

0.01686 

612.7 

12.21 

n/a 

8378 

179.3 

306.4 

12.14 

0.02046 

671.1 

13.49 

n/a 

8425 

188.7 

335.6 

13.42 

0.0223 

991.6 

16.63 

n/a 

8.959 

289.8 

495.8 

16.51 

0.0506 

997 

16.73 

n/a 

8.963 

290.9 

498.5 

16.615 

0.05065 

1150 

19.9 

n/a 

9.072 

325.4 

575 

19.78 

0.05284 

1247 

22.05 

n/a 

9.143 

350.5 

623  5 

21.93 

0.05472 

1292 

23.07 

n/a 

9.176 

367.7 

646 

22.95 

0.5569 

1292 

23.07 

n/a 

9.176 

362.8 

646 

22.95 

0.05569 

1396 

25.51 

n/a 

9.2532 

393.1 

698 

25.38 

0.05795 

1557 

29.44 

n/a 

9.372 

445.9 

778.5 

29.3 

0.06134 

1567 

29  68 

n/a 

9.38 

450 

783.5 

29.55 

0.06212 

1871 

37.58 

n/a 

9.603 

571.2 

935.5 

37.44 

0.0676 

1887 

38.01 

n/a 

9.614 

579.7 

943.5 

37.87 

0.06785 

1900 

38.36 

n/a 

9.624 

584.9 

950 

38.22 

0  06811 

1908 

38.58 

n/a 

9.63 

588.6 

954 

38.43 

0.06824 

2500 

55.51 

n/a 

10.05 

912.8 

1250 

55.33 

0.08296 

3000 

71.19 

n/a 

10.38 

1257 

1500 

71 

0  09485 

3500 

88.12 

n/a 

10.7 

1663 

1750 

87.91 

0.106 

4000 

106.3 

n/a 

11.01 

2132 

2000 

106 

0.1165 

Figure  C.2: 


Spreadsheets  Generated  Using  Iron  Impacts 
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